star-mc

Parallel estimation of Monte Carlo integrators
git clone git://git.meso-star.fr/star-mc.git
Log | Files | Refs | README | LICENSE

test_smc_errors.c (2984B)


      1 /* Copyright (C) 2015-2018, 2021-2023 |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 #include "smc.h"
     17 #include "test_smc_utils.h"
     18 
     19 #include <star/ssp.h>
     20 
     21 static res_T
     22 compute_1
     23   (void* value,
     24    struct ssp_rng* rng,
     25    const unsigned ithread,
     26    const uint64_t irealisation,
     27    void* context)
     28 {
     29   /* pretend that 10% of the realizations fail */
     30   int ok = ssp_rng_canonical(rng) > 0.1;
     31   (void)context, (void)ithread, (void)irealisation;
     32   if(ok) {
     33     SMC_DOUBLE(value) = 1; /* weight is 1 */
     34   } else {
     35     SMC_DOUBLE(value) = 0; /* clear weight */
     36   }
     37   /* a realization has been carried out only if ok */
     38   return ok ? RES_OK : RES_BAD_ARG;
     39 }
     40 
     41 int
     42 main(int argc, char** argv)
     43 {
     44   struct mem_allocator allocator;
     45   struct smc_device* dev;
     46   struct smc_device_create_args args = SMC_DEVICE_CREATE_ARGS_DEFAULT;
     47   struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
     48   struct smc_estimator* estimator;
     49   struct smc_estimator_status status;
     50   (void)argc, (void)argv;
     51 
     52   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
     53 
     54   args.allocator = &allocator;
     55   SMC(device_create(&args, &dev));
     56 
     57   integrator.integrand = &compute_1;
     58   integrator.type = &smc_double;
     59   integrator.max_realisations = 10000;
     60 
     61   SMC(solve(dev, &integrator, NULL, &estimator));
     62   SMC(estimator_get_status(estimator, &status));
     63 
     64   /* result must be 1 with std err = 0 */
     65   CHK(SMC_DOUBLE(status.E) == 1);
     66   CHK(SMC_DOUBLE(status.SE) == 0);
     67 
     68   printf("OK realizations = %lu; KO realizations = %lu\n",
     69 	 (unsigned long)status.N, (unsigned long)status.NF);
     70 
     71   SMC(estimator_ref_put(estimator));
     72 
     73   /* set auto cancel ON */
     74   integrator.max_failures = integrator.max_realisations / 1000;
     75 
     76   SMC(solve(dev, &integrator, NULL, &estimator));
     77   SMC(estimator_get_status(estimator, &status));
     78 
     79   /* it is expected that solve was auto-canceled */
     80   CHK(integrator.max_realisations != status.N);
     81   CHK(status.NF >= integrator.max_failures + 1);
     82 
     83   /* result must be 1 with std err = 0 */
     84   CHK(SMC_DOUBLE(status.E) == 1);
     85   CHK(SMC_DOUBLE(status.SE) == 0);
     86 
     87   printf("OK realizations = %lu; KO realizations = %lu\n",
     88 	 (unsigned long)status.N, (unsigned long)status.NF);
     89 
     90   SMC(device_ref_put(dev));
     91   SMC(estimator_ref_put(estimator));
     92 
     93   check_memory_allocator(&allocator);
     94   mem_shutdown_proxy_allocator(&allocator);
     95   CHK(mem_allocated_size() == 0);
     96   return 0;
     97 }
     98