test_smc_solve.c (4081B)
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 #include <rsys/math.h> 21 22 #include <math.h> 23 24 static res_T 25 rcp_x 26 (void* value, 27 struct ssp_rng* rng, 28 const unsigned ithread, 29 const uint64_t irealisation, 30 void* ctx) 31 { 32 double* result = value; 33 double samp; 34 (void)ithread, (void)irealisation; 35 CHK(value != NULL); 36 CHK(rng != NULL); 37 CHK(ctx == NULL); 38 samp = ssp_rng_uniform_double(rng, 1.0, 4.0); 39 *result = 1.0 / samp * 3; 40 return RES_OK; 41 } 42 43 static res_T 44 cos_x 45 (void* value, 46 struct ssp_rng* rng, 47 const unsigned ithread, 48 const uint64_t irealisation, 49 void* ctx) 50 { 51 float* result = value; 52 double samp; 53 (void)ithread, (void)irealisation; 54 CHK(value != NULL); 55 CHK(rng != NULL); 56 CHK(ctx == (void*)0xC0DE); 57 samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0); 58 *result = (float)(cos(samp) * PI / 2.0); 59 return RES_OK; 60 } 61 62 int 63 main(int argc, char** argv) 64 { 65 struct mem_allocator allocator; 66 struct smc_device* dev; 67 struct smc_device_create_args args = SMC_DEVICE_CREATE_ARGS_DEFAULT; 68 struct smc_integrator integrator = SMC_INTEGRATOR_NULL; 69 struct smc_estimator* estimator; 70 struct smc_estimator_status status; 71 (void)argc, (void)argv; 72 73 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 74 75 args.allocator = &allocator; 76 CHK(smc_device_create(&args, &dev) == RES_OK); 77 78 integrator.integrand = rcp_x; 79 integrator.type = &smc_double; 80 integrator.max_realisations = 10000; 81 82 CHK(smc_solve(NULL, NULL, NULL, NULL) == RES_BAD_ARG); 83 CHK(smc_solve(dev, NULL, NULL, NULL) == RES_BAD_ARG); 84 CHK(smc_solve(NULL, &integrator, NULL, NULL) == RES_BAD_ARG); 85 CHK(smc_solve(dev, &integrator, NULL, NULL) == RES_BAD_ARG); 86 CHK(smc_solve(NULL, NULL, NULL, &estimator) == RES_BAD_ARG); 87 CHK(smc_solve(dev, NULL, NULL, &estimator) == RES_BAD_ARG); 88 CHK(smc_solve(NULL, &integrator, NULL, &estimator) == RES_BAD_ARG); 89 CHK(smc_solve(dev, &integrator, NULL, &estimator) == RES_OK); 90 91 CHK(smc_estimator_get_status(NULL, NULL) == RES_BAD_ARG); 92 CHK(smc_estimator_get_status(estimator, NULL) == RES_BAD_ARG); 93 CHK(smc_estimator_get_status(NULL, &status) == RES_BAD_ARG); 94 CHK(smc_estimator_get_status(estimator, &status) == RES_OK); 95 printf("Integral[1, 4] 1/x dx = %g; E = %g; SE = %g\n", 96 log(4.0) - log(1.0), SMC_DOUBLE(status.E), SMC_DOUBLE(status.SE)); 97 CHK(eq_eps 98 (log(4.0) - log(1.0), SMC_DOUBLE(status.E), 2.0 * SMC_DOUBLE(status.SE))); 99 CHK(smc_estimator_ref_put(estimator) == RES_OK); 100 101 integrator.type = NULL; 102 CHK(smc_solve(dev, &integrator, NULL, &estimator) == RES_BAD_ARG); 103 integrator.type = &smc_double; 104 integrator.integrand = NULL; 105 CHK(smc_solve(dev, &integrator, NULL, &estimator) == RES_BAD_ARG); 106 107 integrator.integrand = cos_x; 108 integrator.type = &smc_float; 109 CHK(smc_solve(dev, &integrator, (void*)0xC0DE, &estimator) == RES_OK); 110 CHK(smc_estimator_get_status(estimator, &status) == RES_OK); 111 printf("Integral[-PI/4, PI/4] cos(x) dx = %f; E = %f; SE = %f\n", 112 2*sin(PI/4.0), SMC_FLOAT(status.E), SMC_FLOAT(status.SE)); 113 CHK(eq_eps 114 ((float)2*sin(PI/4.0), SMC_FLOAT(status.E), 2 * SMC_FLOAT(status.SE))); 115 CHK(smc_device_ref_put(dev) == RES_OK); 116 117 CHK(smc_estimator_ref_put(estimator) == RES_OK); 118 119 check_memory_allocator(&allocator); 120 mem_shutdown_proxy_allocator(&allocator); 121 CHK(mem_allocated_size() == 0); 122 return 0; 123 } 124