commit 86820bac064a09b79faff66693163a579c645672
parent 89334116bca62fd013c5f005528459a1da9db518
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Jul 2015 12:16:42 +0200
Test the smc_integrate function
Diffstat:
8 files changed, 182 insertions(+), 15 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -124,6 +124,7 @@ if(NOT NO_TEST)
new_test(test_smc_device)
new_test(test_smc_solve ${MATH_LIB})
+ new_test(test_smc_integrate ${MATH_LIB})
if(NOT TARGET Star3D)
rcmake_copy_runtime_libraries(test_smc_solve)
diff --git a/src/smc.h b/src/smc.h
@@ -104,13 +104,16 @@ struct smc_integrator {
size_t max_steps; /* Maximum # of experiments */
};
+static const struct smc_integrator SMC_INTEGRATOR_NULL =
+{ NULL, NULL, SIZE_MAX };
+
struct smc_integrator_raw {
void (*integrand)(void* accum, struct ssp_rng* rng, void* ctx);
const struct smc_type_raw* type;
size_t max_steps; /* Maximum # of accumulation */
};
-static const struct smc_integrator SMC_INTEGRATOR_NULL =
+static const struct smc_integrator_raw SMC_INTEGRATOR_RAW_NULL =
{ NULL, NULL, SIZE_MAX };
/* Forward declaration of opaque types */
diff --git a/src/smc_accumulator.c b/src/smc_accumulator.c
@@ -162,7 +162,7 @@ smc_integrate
#pragma omp parallel for schedule(static)
for(i=0; i < (int64_t)integrator->max_steps; ++i) {
const int ithread = omp_get_thread_num();
- integrator->type->clear(vals[i]);
+ integrator->type->clear(vals[ithread]);
integrator->integrand(vals[ithread], dev->rngs[ithread], ctx);
integrator->type->add(accums[ithread], accums[ithread], vals[ithread]);
++nsamples[ithread];
diff --git a/src/test_smc_device.c b/src/test_smc_device.c
@@ -1,8 +1,5 @@
/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
*
- * This software is a computer program whose purpose is to manage the
- * statistical estimation of a function.
- *
* This software is governed by the CeCILL license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/or redistribute the software under the terms of the CeCILL
diff --git a/src/test_smc_integrate.c b/src/test_smc_integrate.c
@@ -0,0 +1,175 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include "smc.h"
+#include "test_smc_utils.h"
+
+#include <star/ssp.h>
+
+#include <math.h>
+
+static void*
+dbl2_create(struct mem_allocator* allocator, void* ctx)
+{
+ NCHECK(allocator, NULL);
+ CHECK(ctx, (void*)0xC0DE);
+ return MEM_CALLOC(allocator, 2, sizeof(double));
+}
+
+static void
+dbl2_destroy(struct mem_allocator* allocator, void* data)
+{
+ NCHECK(allocator, NULL);
+ NCHECK(data, NULL);
+ MEM_RM(allocator, data);
+}
+
+static void
+dbl2_clear(void* data)
+{
+ ((double*)data)[0] = ((double*)data)[1] = 0.0;
+}
+
+static void
+dbl2_add(void* result, const void* op0, const void* op1)
+{
+ int i;
+ FOR_EACH(i, 0, 2) {
+ ((double*)result)[i] = ((double*)op0)[i] + ((double*)op1)[i];
+ }
+}
+
+static const struct smc_type_raw smc_accum_dbl2 = {
+ dbl2_create,
+ dbl2_destroy,
+ dbl2_clear,
+ dbl2_add
+};
+
+static void
+rcp_x(void* value, struct ssp_rng* rng, void* ctx)
+{
+ double* result = value;
+ double samp;
+ NCHECK(value, NULL);
+ NCHECK(rng, NULL);
+ CHECK(ctx, (void*)0xC0DE);
+ samp = ssp_rng_uniform_double(rng, 1.0, 4.0);
+ result[0] = 1.0 / samp * 3;
+ result[1] = result[0]*result[0];
+}
+
+static void
+cos_x(void* value, struct ssp_rng* rng, void* ctx)
+{
+ double* result = value;
+ double samp;
+ NCHECK(value, NULL);
+ NCHECK(rng, NULL);
+ CHECK(ctx, (void*)0xC0DE);
+ samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0);
+ result[0] = cos(samp) * PI / 2.0;
+ result[1] = result[0]*result[0];
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct smc_device* dev;
+ struct smc_integrator_raw integrator = SMC_INTEGRATOR_RAW_NULL;
+ struct smc_accumulator* accum;
+ struct smc_accumulator_status status;
+ double* val;
+ double E, V, SE;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+ CHECK(smc_device_create(NULL, &allocator, SMC_NTHREADS_DEFAULT, &dev), RES_OK);
+
+ integrator.integrand = rcp_x;
+ integrator.type = &smc_accum_dbl2;
+ integrator.max_steps = 10000;
+
+ CHECK(smc_integrate(NULL, NULL, (void*)0xC0DE, NULL), RES_BAD_ARG);
+ CHECK(smc_integrate(dev, NULL, (void*)0xC0DE, NULL), RES_BAD_ARG);
+ CHECK(smc_integrate(NULL, &integrator, (void*)0xC0DE, NULL), RES_BAD_ARG);
+ CHECK(smc_integrate(dev, &integrator, (void*)0xC0DE, NULL), RES_BAD_ARG);
+ CHECK(smc_integrate(NULL, NULL, (void*)0xC0DE, &accum), RES_BAD_ARG);
+ CHECK(smc_integrate(dev, NULL, (void*)0xC0DE, &accum), RES_BAD_ARG);
+ CHECK(smc_integrate(NULL, &integrator, (void*)0xC0DE, &accum), RES_BAD_ARG);
+ CHECK(smc_integrate(dev, &integrator, (void*)0xC0DE, &accum), RES_OK);
+
+ CHECK(smc_accumulator_get_status(NULL, NULL), RES_BAD_ARG);
+ CHECK(smc_accumulator_get_status(accum, NULL), RES_BAD_ARG);
+ CHECK(smc_accumulator_get_status(NULL, &status), RES_BAD_ARG);
+ CHECK(smc_accumulator_get_status(accum, &status), RES_OK);
+
+ val = status.value;
+ E = val[0] / (double)status.N;
+ V = val[1] / (double)status.N - (val[0]*val[0]) / (double)(status.N * status.N);
+ SE = sqrt(V / (double)status.N);
+
+ printf("Integral[1, 4] 1/x dx = %g; E = %g; SE = %g\n",
+ log(4.0) - log(1.0), E, SE);
+ CHECK(eq_eps(log(4.0) - log(1.0), E, 2.0 * SE), 1);
+
+ CHECK(smc_accumulator_ref_get(NULL), RES_BAD_ARG);
+ CHECK(smc_accumulator_ref_get(accum), RES_OK);
+ CHECK(smc_accumulator_ref_put(NULL), RES_BAD_ARG);
+ CHECK(smc_accumulator_ref_put(accum), RES_OK);
+ CHECK(smc_accumulator_ref_put(accum), RES_OK);
+
+ integrator.type = NULL;
+ CHECK(smc_integrate(dev, &integrator, NULL, &accum), RES_BAD_ARG);
+ integrator.type = &smc_accum_dbl2;
+ integrator.integrand = NULL;
+ CHECK(smc_integrate(dev, &integrator, NULL, &accum), RES_BAD_ARG);
+
+ integrator.integrand = cos_x;
+ CHECK(smc_integrate(dev, &integrator, (void*)0xC0DE, &accum), RES_OK);
+ CHECK(smc_accumulator_get_status(accum, &status), RES_OK);
+
+ val = status.value;
+ E = val[0] / (double)status.N;
+ V = val[1] / (double)status.N - (val[0]*val[0]) / (double)(status.N * status.N);
+ SE = sqrt(V / (double)status.N);
+
+ printf("Integral[-PI/4, PI/4] cos(x) dx = %g; E = %g; SE = %g\n",
+ 2*sin(PI/4.0), E, SE);
+ CHECK(eq_eps(2*sin(PI/4.0), E, 2 * SE), 1);
+
+ CHECK(smc_accumulator_ref_put(accum), RES_OK);
+ CHECK(smc_device_ref_put(dev), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
+
diff --git a/src/test_smc_light_path.c b/src/test_smc_light_path.c
@@ -1,8 +1,5 @@
/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
*
- * This software is a computer program whose purpose is to manage the
- * statistical estimation of a function.
- *
* This software is governed by the CeCILL license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/or redistribute the software under the terms of the CeCILL
diff --git a/src/test_smc_solve.c b/src/test_smc_solve.c
@@ -1,8 +1,5 @@
/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
*
- * This software is a computer program whose purpose is to manage the
- * statistical estimation of a function.
- *
* This software is governed by the CeCILL license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/or redistribute the software under the terms of the CeCILL
@@ -92,7 +89,7 @@ main(int argc, char** argv)
CHECK(smc_estimator_get_status(estimator, NULL), RES_BAD_ARG);
CHECK(smc_estimator_get_status(NULL, &status), RES_BAD_ARG);
CHECK(smc_estimator_get_status(estimator, &status), RES_OK);
- printf("Integral[1, 4] 1/x dx = %f; E = %f; SE = %f\n",
+ printf("Integral[1, 4] 1/x dx = %g; E = %g; SE = %g\n",
log(4.0) - log(1.0), SMC_DOUBLE(status.E), SMC_DOUBLE(status.SE));
CHECK(eq_eps
(log(4.0) - log(1.0), SMC_DOUBLE(status.E), 2.0 * SMC_DOUBLE(status.SE)), 1);
diff --git a/src/test_smc_utils.h b/src/test_smc_utils.h
@@ -1,8 +1,5 @@
/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
*
- * This software is a computer program whose purpose is to manage the
- * statistical estimation of a function.
- *
* This software is governed by the CeCILL license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/or redistribute the software under the terms of the CeCILL
