test_smc_doubleN.c (3544B)
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 #include <math.h> 22 23 enum { 24 WEIGHT, 25 WEIGHT_SENSIBILITY_ALPHA, 26 WEIGHTS_COUNT 27 }; 28 29 struct alpha { double value; }; 30 31 static res_T 32 cos_x 33 (void* value, 34 struct ssp_rng* rng, 35 const unsigned ithread, 36 const uint64_t irealisation, 37 void* context) 38 { 39 double* result = SMC_DOUBLEN(value); 40 double samp; 41 struct smc_doubleN_context* ctx = context; 42 const struct alpha* alpha = ctx->integrand_data; 43 (void)ithread, (void)irealisation; 44 CHK(value != NULL); 45 CHK(rng != NULL); 46 samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0); 47 48 result[WEIGHT] = cos(alpha->value*samp) * PI / 2.0; 49 result[WEIGHT_SENSIBILITY_ALPHA] = -samp * sin(alpha->value*samp) * PI/2.0; 50 51 return RES_OK; 52 } 53 54 int 55 main(int argc, char** argv) 56 { 57 struct mem_allocator allocator; 58 struct smc_device* dev; 59 struct smc_device_create_args args = SMC_DEVICE_CREATE_ARGS_DEFAULT; 60 struct smc_integrator integrator = SMC_INTEGRATOR_NULL; 61 struct smc_estimator* estimator; 62 struct smc_estimator_status status; 63 struct smc_doubleN_context ctx; 64 struct alpha alpha; 65 const double a = 0.4; 66 double result; 67 (void)argc, (void)argv; 68 69 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 70 71 args.allocator = &allocator; 72 CHK(smc_device_create(&args, &dev) == RES_OK); 73 74 integrator.integrand = cos_x; 75 integrator.type = &smc_doubleN; 76 integrator.max_realisations = 100000; 77 alpha.value = a; 78 ctx.count = WEIGHTS_COUNT; 79 ctx.integrand_data = α 80 81 CHK(smc_solve(dev, &integrator, NULL, &estimator) == RES_MEM_ERR); 82 ctx.count = 0; 83 CHK(smc_solve(dev, &integrator, &ctx, &estimator) == RES_MEM_ERR); 84 ctx.count = WEIGHTS_COUNT; 85 CHK(smc_solve(dev, &integrator, &ctx, &estimator) == RES_OK); 86 CHK(smc_estimator_get_status(estimator, &status) == RES_OK); 87 88 /* Check the estimated result */ 89 result = 2*sin(a*PI/4.0) / a; 90 printf("Integral[-PI/4, PI/4] cos(%g*x) dx = %g ~ %g +/- %g\n", 91 a, result, 92 SMC_DOUBLEN(status.E)[WEIGHT], 93 SMC_DOUBLEN(status.SE)[WEIGHT]); 94 CHK(eq_eps 95 (result, 96 SMC_DOUBLEN(status.E)[WEIGHT], 97 SMC_DOUBLEN(status.SE)[WEIGHT])); 98 99 /* Check the estimated sensibility in alpha */ 100 result = 2*(a*PI/4.0*cos(a*PI/4.0) - sin(a*PI/4.0)) / (a*a); 101 printf("Integral'[-PI/4, PI/4] cos(%g*x) dx = %g ~ %g +/- %g\n", 102 a, result, 103 SMC_DOUBLEN(status.E)[WEIGHT_SENSIBILITY_ALPHA], 104 SMC_DOUBLEN(status.SE)[WEIGHT_SENSIBILITY_ALPHA]); 105 CHK(eq_eps 106 (result, 107 SMC_DOUBLEN(status.E)[WEIGHT_SENSIBILITY_ALPHA], 108 SMC_DOUBLEN(status.SE)[WEIGHT_SENSIBILITY_ALPHA])); 109 110 CHK(smc_estimator_ref_put(estimator) == RES_OK); 111 CHK(smc_device_ref_put(dev) == RES_OK); 112 113 check_memory_allocator(&allocator); 114 mem_shutdown_proxy_allocator(&allocator); 115 CHK(mem_allocated_size() == 0); 116 return 0; 117 }