commit 89334116bca62fd013c5f005528459a1da9db518
parent 81b041b45a2fa0bdfe2da2cbceb9c4ffb1c0e76c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Jul 2015 11:42:25 +0200
Add the smc_integrate function
This function performs a raw integration, i.e. the estimation of the
integration must be performed by the caller with respect to accumulated
values provided by the status of the smc_accumulator.
Diffstat:
4 files changed, 285 insertions(+), 2 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -64,6 +64,7 @@ set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SMC_FILES_SRC
+ smc_accumulator.c
smc_device.c
smc_integrator.c
smc_type.c)
diff --git a/src/smc.h b/src/smc.h
@@ -74,10 +74,18 @@ struct smc_type {
void (*sqrt)(void* result, const void* value);
};
+/* Type for raw integration, i.e. whithout estimator. Experiment results are
+ * simply accumulated */
+struct smc_type_raw {
+ void* (*create)(struct mem_allocator* allocator, void* ctx);
+ void (*destroy)(struct mem_allocator*, void* accum);
+ void (*clear)(void* accum);
+ void (*add)(void* result, const void* op0, const void* op1);
+};
+
static const struct smc_type SMC_TYPE_NULL =
{ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
-
struct smc_estimator_status {
void* E; /* Expected value */
void* V; /* Variance */
@@ -85,18 +93,30 @@ struct smc_estimator_status {
size_t N; /* Samples count */
};
+struct smc_accumulator_status {
+ void* value; /* Accumulated value */
+ size_t N; /* Samples count */
+};
+
struct smc_integrator {
void (*integrand)(void* value, struct ssp_rng* rng, void* ctx);
const struct smc_type* type;
size_t max_steps; /* Maximum # of experiments */
};
+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 =
{ NULL, NULL, SIZE_MAX };
/* Forward declaration of opaque types */
struct smc_device; /* Entry point of the library */
struct smc_estimator; /* Estimator of an integrator */
+struct smc_accumulator; /* Accumulated values of a integration */
BEGIN_DECLS
@@ -109,7 +129,7 @@ SMC_API const struct smc_type smc_double;
#define SMC_DOUBLE(Val) (*(double*)(Val))
/*******************************************************************************
- * API deinition
+ * Device API
******************************************************************************/
SMC_API res_T
smc_device_create
@@ -126,6 +146,9 @@ SMC_API res_T
smc_device_ref_put
(struct smc_device* dev);
+/*******************************************************************************
+ * Integration API
+ ******************************************************************************/
SMC_API res_T
smc_solve
(struct smc_device* dev,
@@ -155,6 +178,29 @@ smc_estimator_get_status
(struct smc_estimator* estimator,
struct smc_estimator_status* status);
+/*******************************************************************************
+ * Accumulation API
+ ******************************************************************************/
+SMC_API res_T
+smc_integrate
+ (struct smc_device* dev,
+ struct smc_integrator_raw* integrator,
+ void* ctx,
+ struct smc_accumulator** accumulator);
+
+SMC_API res_T
+smc_accumulator_ref_get
+ (struct smc_accumulator* accumulator);
+
+SMC_API res_T
+smc_accumulator_ref_put
+ (struct smc_accumulator* accumulator);
+
+SMC_API res_T
+smc_accumulator_get_status
+ (struct smc_accumulator* accumulator,
+ struct smc_accumulator_status* status);
+
END_DECLS
#endif /* SMC_H */
diff --git a/src/smc_accumulator.c b/src/smc_accumulator.c
@@ -0,0 +1,225 @@
+/* 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 "smc_device_c.h"
+#include "smc_type_c.h"
+
+#include <rsys/mem_allocator.h>
+
+#include <omp.h>
+
+struct smc_accumulator {
+ struct smc_type_raw type;
+ struct smc_accumulator_status status;
+ struct smc_device* dev;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+accumulator_create
+ (struct smc_device* dev,
+ const struct smc_type_raw* type,
+ void* ctx,
+ struct smc_accumulator** out_accum)
+{
+ struct smc_accumulator* accum = NULL;
+ res_T res = RES_OK;
+ ASSERT(out_accum && dev && type);
+
+ accum = MEM_CALLOC(dev->allocator, 1, sizeof(struct smc_accumulator));
+ if(!accum) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ SMC(device_ref_get(dev));
+ accum->dev = dev;
+ ref_init(&accum->ref);
+
+ accum->status.value = type->create(dev->allocator, ctx);
+ if(!accum->status.value) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ type->clear(accum->status.value);
+
+ accum->status.N = 0;
+ accum->type = *type;
+
+exit:
+ *out_accum = accum;
+ return res;
+error:
+ if(accum) {
+ SMC(accumulator_ref_put(accum));
+ accum = NULL;
+ }
+ goto exit;
+}
+
+static char
+check_integrator_raw(struct smc_integrator_raw* integrator)
+{
+ ASSERT(integrator);
+ return integrator->integrand
+ && integrator->type
+ && integrator->max_steps
+ && check_type_raw(integrator->type);
+}
+
+static void
+accumulator_release(ref_T* ref)
+{
+ struct smc_accumulator* accum;
+ struct smc_device* dev;
+ ASSERT(ref);
+
+ accum = CONTAINER_OF(ref, struct smc_accumulator, ref);
+ dev = accum->dev;
+ if(accum->status.value)
+ accum->type.destroy(dev->allocator, accum->status.value);
+ MEM_RM(dev->allocator, accum);
+ SMC(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+smc_integrate
+ (struct smc_device* dev,
+ struct smc_integrator_raw* integrator,
+ void* ctx,
+ struct smc_accumulator** out_accum)
+{
+ struct smc_accumulator* accum = NULL;
+ int64_t i;
+ size_t nthreads = 0;
+ void** vals = NULL;
+ void** accums = NULL;
+ size_t* nsamples = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !integrator || !out_accum || !check_integrator_raw(integrator)) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ nthreads = device_get_threads_count(dev);
+
+ /* Create and initializdd per thread memory variables */
+ nsamples = MEM_CALLOC(dev->allocator, nthreads, sizeof(size_t));
+ if(!nsamples) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ vals = MEM_CALLOC(dev->allocator, nthreads, sizeof(void*));
+ if(!vals) { res = RES_MEM_ERR; goto error; }
+ accums = MEM_CALLOC(dev->allocator, nthreads, sizeof(void*));
+ FOR_EACH(i, 0, (int64_t)nthreads) {
+ vals[i] = integrator->type->create(dev->allocator, ctx);
+ if(!vals[i]) { res = RES_MEM_ERR; goto error; }
+ accums[i] = integrator->type->create(dev->allocator, ctx);
+ if(!accums[i]) { res = RES_MEM_ERR; goto error; }
+ integrator->type->clear(accums[i]);
+ }
+
+ /* Create the accumulator */
+ res = accumulator_create(dev, integrator->type, ctx, &accum);
+ if(res != RES_OK) goto error;
+
+ /* Parallel evaluation of the integrator */
+ #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->integrand(vals[ithread], dev->rngs[ithread], ctx);
+ integrator->type->add(accums[ithread], accums[ithread], vals[ithread]);
+ ++nsamples[ithread];
+ }
+
+ /* Merge the parallel estimation into the final estimator */
+ FOR_EACH(i, 0, (int64_t)nthreads) {
+ accum->status.N += nsamples[i];
+ integrator->type->add(accum->status.value, accum->status.value, accums[i]);
+ }
+
+exit:
+ if(nsamples) MEM_RM(dev->allocator, nsamples);
+ if(vals) {
+ FOR_EACH(i, 0, (int64_t)nthreads)
+ if(vals[i]) integrator->type->destroy(dev->allocator, vals[i]);
+ MEM_RM(dev->allocator, vals);
+ }
+ if(accums) {
+ FOR_EACH(i, 0, (int64_t)nthreads)
+ if(accums[i]) integrator->type->destroy(dev->allocator, accums[i]);
+ MEM_RM(dev->allocator, accums);
+ }
+ if(out_accum) *out_accum = accum;
+ return res;
+
+error:
+ if(accum) {
+ SMC(accumulator_ref_put(accum));
+ accum = NULL;
+ }
+ goto exit;
+}
+
+res_T
+smc_accumulator_ref_get(struct smc_accumulator* accum)
+{
+ if(!accum) return RES_BAD_ARG;
+ ref_get(&accum->ref);
+ return RES_OK;
+}
+
+res_T
+smc_accumulator_ref_put(struct smc_accumulator* accum)
+{
+ if(!accum) return RES_BAD_ARG;
+ ref_put(&accum->ref, accumulator_release);
+ return RES_OK;
+}
+
+res_T
+smc_accumulator_get_status
+ (struct smc_accumulator* accum,
+ struct smc_accumulator_status* status)
+{
+ if(!accum || !status) return RES_BAD_ARG;
+ *status = accum->status;
+ return RES_OK;
+}
+
diff --git a/src/smc_type_c.h b/src/smc_type_c.h
@@ -37,6 +37,7 @@
static INLINE char
check_type(const struct smc_type* type)
{
+ ASSERT(type);
return type->create != NULL
&& type->destroy != NULL
&& type->set != NULL
@@ -48,5 +49,15 @@ check_type(const struct smc_type* type)
&& type->sqrt != NULL;
}
+static INLINE char
+check_type_raw(const struct smc_type_raw* type)
+{
+ ASSERT(type);
+ return type->create != NULL
+ && type->destroy != NULL
+ && type->clear != NULL
+ && type->add != NULL;
+}
+
#endif /* SMC_TYPE_C_H */
