commit 7e3cfe4cf5e8e5f4a6f0c56f96b7f08a7c1010dd
parent 41003e0912f6baf80ab1d32bb0a968c0891990d6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 7 Apr 2016 16:08:55 +0200
Add and test the smc_doubleN builtin type
Diffstat:
4 files changed, 355 insertions(+), 2 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -66,6 +66,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SMC_FILES_SRC
smc_accumulator.c
smc_device.c
+ smc_doubleN.c
smc_estimator.c
smc_type.c)
set(SMC_FILES_INC_API smc.h)
@@ -89,7 +90,6 @@ set_target_properties(smc PROPERTIES
if(OPENMP_FOUND)
set_target_properties(smc PROPERTIES COMPILE_FLAGS ${OpenMP_C_FLAGS})
-
if(CMAKE_COMPILER_IS_GNUCC)
set_target_properties(smc PROPERTIES LINK_FLAGS ${OpenMP_C_FLAGS})
endif()
@@ -124,6 +124,7 @@ if(NOT NO_TEST)
new_test(test_smc_device)
new_test(test_smc_solve ${MATH_LIB})
+ new_test(test_smc_doubleN ${MATH_LIB})
new_test(test_smc_integrate ${MATH_LIB})
if(NOT TARGET Star3D)
diff --git a/src/smc.h b/src/smc.h
@@ -127,10 +127,20 @@ BEGIN_DECLS
/* Pre-declared SMC types */
SMC_API const struct smc_type smc_float;
SMC_API const struct smc_type smc_double;
+SMC_API const struct smc_type smc_doubleN;
+
+/* Context of a list of double precision floating point data. Must be set as
+ * the context parameter of the smc_solve function when the integrator type is
+ * smc_doubleN. */
+struct smc_doubleN_context {
+ size_t count; /* Number of floating point values */
+ void* integrand_data; /* User defined data used by the integrand */
+};
-/* Syntactic sugar macros */
+/* Syntactic sugar macros and functions */
#define SMC_FLOAT(Val) (*(float*)(Val))
#define SMC_DOUBLE(Val) (*(double*)(Val))
+SMC_API double* SMC_DOUBLEN(void* val);
/*******************************************************************************
* Device API
diff --git a/src/smc_doubleN.c b/src/smc_doubleN.c
@@ -0,0 +1,216 @@
+/* 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
+ * 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 <rsys/dynamic_array_double.h>
+#include <math.h>
+
+/*******************************************************************************
+ * smc_doubleN functions
+ ******************************************************************************/
+static void
+smc_doubleN_destroy(struct mem_allocator* allocator, void* data)
+{
+ struct darray_double* dblN = data;
+ ASSERT(data);
+ darray_double_release(dblN);
+ MEM_RM(allocator, dblN);
+}
+
+static void*
+smc_doubleN_create(struct mem_allocator* allocator, void* context)
+{
+ struct smc_doubleN_context* ctx = context;
+ struct darray_double* dblN = NULL;
+ res_T res = RES_OK;
+ ASSERT(allocator);
+
+ if(!ctx || !ctx->count) goto error;
+
+ dblN = MEM_ALLOC(allocator, sizeof(struct darray_double));
+ if(!dblN) goto error;
+
+ darray_double_init(allocator, dblN);
+ res = darray_double_resize(dblN, ctx->count);
+ if(res != RES_OK) goto error;
+
+exit:
+ return dblN;
+error:
+ if(dblN) {
+ smc_doubleN_destroy(allocator, dblN);
+ dblN = NULL;
+ }
+ goto exit;
+}
+
+static void
+smc_doubleN_set(void* result, const void* value)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* src = value;
+ size_t i, n;
+ ASSERT(dst && src);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(src)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] = darray_double_cdata_get(src)[i];
+ }
+}
+
+static void
+smc_doubleN_zero(void* result)
+{
+ struct darray_double* dblN = result;
+ ASSERT(result);
+ memset(darray_double_data_get(dblN), 0,
+ darray_double_size_get(dblN)*sizeof(double));
+}
+
+static void
+smc_doubleN_add(void* result, const void* op0, const void* op1)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* a = op0;
+ const struct darray_double* b = op1;
+ size_t i, n;
+ ASSERT(result && op0 && op1);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(a) || n != darray_double_size_get(b)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] =
+ darray_double_cdata_get(a)[i] + darray_double_cdata_get(b)[i];
+ }
+}
+
+static void
+smc_doubleN_sub(void* result, const void* op0, const void* op1)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* a = op0;
+ const struct darray_double* b = op1;
+ size_t i, n;
+ ASSERT(result && op0 && op1);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(a) || n != darray_double_size_get(b)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] =
+ darray_double_cdata_get(a)[i] - darray_double_cdata_get(b)[i];
+ }
+}
+
+static void
+smc_doubleN_mul(void* result, const void* op0, const void* op1)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* a = op0;
+ const struct darray_double* b = op1;
+ size_t i, n;
+ ASSERT(result && op0 && op1);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(a) || n != darray_double_size_get(b)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] =
+ darray_double_cdata_get(a)[i] * darray_double_cdata_get(b)[i];
+ }
+}
+
+static void
+smc_doubleN_divi(void* result, const void* op0, const size_t op1)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* a = op0;
+ size_t i, n;
+ ASSERT(result && op0 && op1);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(a)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] = darray_double_cdata_get(a)[i]/(double)op1;
+ }
+}
+
+static void
+smc_doubleN_sqrt(void* result, const void* value)
+{
+ struct darray_double* dst = result;
+ const struct darray_double* src = value;
+ size_t i, n;
+ ASSERT(result && value);
+
+ n = darray_double_size_get(dst);
+ if(n != darray_double_size_get(src)) {
+ FATAL("The vectors must have the same dimension.\n");
+ }
+ FOR_EACH(i, 0, n) {
+ darray_double_data_get(dst)[i] = sqrt(darray_double_cdata_get(src)[i]);
+ }
+}
+
+/*******************************************************************************
+ * Exported helper function
+ ******************************************************************************/
+double*
+SMC_DOUBLEN(void* val)
+{
+ ASSERT(val);
+ return darray_double_data_get(val);
+}
+
+/*******************************************************************************
+ * Exported type
+ ******************************************************************************/
+const struct smc_type smc_doubleN = {
+ smc_doubleN_create,
+ smc_doubleN_destroy,
+ smc_doubleN_set,
+ smc_doubleN_zero,
+ smc_doubleN_add,
+ smc_doubleN_sub,
+ smc_doubleN_mul,
+ smc_doubleN_divi,
+ smc_doubleN_sqrt
+};
+
diff --git a/src/test_smc_doubleN.c b/src/test_smc_doubleN.c
@@ -0,0 +1,126 @@
+/* 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
+ * 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>
+
+enum {
+ WEIGHT,
+ WEIGHT_SENSIBILITY_ALPHA,
+ WEIGHTS_COUNT
+};
+
+struct alpha { double value; };
+
+static void
+cos_x(void* value, struct ssp_rng* rng, void* context)
+{
+ double* result = SMC_DOUBLEN(value);
+ double samp;
+ struct smc_doubleN_context* ctx = context;
+ const struct alpha* alpha = ctx->integrand_data;
+ NCHECK(value, NULL);
+ NCHECK(rng, NULL);
+ samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0);
+
+ result[WEIGHT] =
+ cos(alpha->value*samp) * PI / 2.0;
+ result[WEIGHT_SENSIBILITY_ALPHA] =
+ -samp * sin(alpha->value*samp) * PI/2.0;
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct smc_device* dev;
+ struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
+ struct smc_estimator* estimator;
+ struct smc_estimator_status status;
+ struct smc_doubleN_context ctx;
+ struct alpha alpha;
+ const double a = 0.4;
+ double result;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(smc_device_create(NULL, &allocator, SMC_NTHREADS_DEFAULT, NULL, &dev),
+ RES_OK);
+
+ integrator.integrand = cos_x;
+ integrator.type = &smc_doubleN;
+ integrator.max_steps = 100000;
+ alpha.value = a;
+ ctx.count = WEIGHTS_COUNT;
+ ctx.integrand_data = α
+
+ CHECK(smc_solve(dev, &integrator, NULL, &estimator), RES_MEM_ERR);
+ ctx.count = 0;
+ CHECK(smc_solve(dev, &integrator, &ctx, &estimator), RES_MEM_ERR);
+ ctx.count = WEIGHTS_COUNT;
+ CHECK(smc_solve(dev, &integrator, &ctx, &estimator), RES_OK);
+ CHECK(smc_estimator_get_status(estimator, &status), RES_OK);
+
+ /* Check the estimated result */
+ result = 2*sin(a*PI/4.0) / a;
+ printf("Integral[-PI/4, PI/4] cos(%g*x) dx = %g ~ %g +/- %g\n",
+ a, result,
+ SMC_DOUBLEN(status.E)[WEIGHT],
+ SMC_DOUBLEN(status.SE)[WEIGHT]);
+ CHECK(eq_eps
+ (result,
+ SMC_DOUBLEN(status.E)[WEIGHT],
+ SMC_DOUBLEN(status.SE)[WEIGHT]), 1);
+
+ /* Check the estimated sensibility in alpha */
+ result = 2*(a*PI/4.0*cos(a*PI/4.0) - sin(a*PI/4.0)) / (a*a);
+ printf("Integral'[-PI/4, PI/4] cos(%g*x) dx = %g ~ %g +/- %g\n",
+ a, result,
+ SMC_DOUBLEN(status.E)[WEIGHT_SENSIBILITY_ALPHA],
+ SMC_DOUBLEN(status.SE)[WEIGHT_SENSIBILITY_ALPHA]);
+ CHECK(eq_eps
+ (result,
+ SMC_DOUBLEN(status.E)[WEIGHT_SENSIBILITY_ALPHA],
+ SMC_DOUBLEN(status.SE)[WEIGHT_SENSIBILITY_ALPHA]), 1);
+
+ CHECK(smc_estimator_ref_put(estimator), 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;
+}
