star-mc

Parallel estimation of Monte Carlo integrators
git clone git://git.meso-star.fr/star-mc.git
Log | Files | Refs | README | LICENSE

commit 818b520e038f667730618f02cae7a6c9b0247ebc
parent 2bfea76647655ea5db95fa080ee7f35569209236
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  1 Apr 2015 11:39:39 +0200

Implement the smc_estimator_get_status function

Diffstat:
Mcmake/CMakeLists.txt | 3+++
Msrc/smc.h | 20+++++++++++++++++++-
Msrc/smc_integrator.c | 97++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------
Msrc/smc_type.c | 11++++++++++-
Msrc/smc_type_c.h | 3++-
Msrc/test_smc_device.c | 1-
6 files changed, 111 insertions(+), 24 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -68,6 +68,9 @@ set_target_properties(smc PROPERTIES SOVERSION ${VERSION_MAJOR}) target_link_libraries(smc RSys) +if(CMAKE_COMPILER_IS_GNUCC) + target_link_libraries(smc m) +endif() rcmake_setup_devel(smc StarMC ${VERSION} star/smc_version.h) diff --git a/src/smc.h b/src/smc.h @@ -67,8 +67,19 @@ struct smc_type { void (*sub)(void* result, const void* op0, const void* op1); void (*mul)(void* result, const void* op0, const void* op1); void (*divi)(void* result, const void* op0, const unsigned long op1); + void (*sqrt)(void* result, const void* value); }; +struct smc_estimator_status { + void* expected_value; + void* variance; + void* standard_error; + unsigned long samples_count; +}; + +static const struct smc_type smc_type_null = +{ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL }; + /* Forward declaration of opaque types */ struct smc_device; /* Entry point of the library */ struct smc_estimator; /* Estimator of an integrator */ @@ -76,7 +87,9 @@ struct smc_estimator; /* Estimator of an integrator */ BEGIN_DECLS /* Pre-declared SMC types */ -extern const struct smc_type smc_float; +SMC_API const struct smc_type smc_float; + +#define SMC_FLOAT(Val) (*(const float*)(Val)) /******************************************************************************* * API deinition @@ -111,6 +124,11 @@ SMC_API res_T smc_estimator_ref_put (struct smc_estimator* estimator); +SMC_API res_T +smc_estimator_get_status + (struct smc_estimator* estimator, + struct smc_estimator_status* status); + END_DECLS #endif /* SMC_H */ diff --git a/src/smc_integrator.c b/src/smc_integrator.c @@ -38,10 +38,13 @@ #include <limits.h> struct smc_estimator { + struct smc_type type; void* value; void* square_value; unsigned long nsamples; + struct smc_estimator_status status; + struct smc_device* dev; ref_T ref; }; @@ -58,8 +61,16 @@ estimator_release(ref_T* ref) estimator = CONTAINER_OF(ref, struct smc_estimator, ref); dev = estimator->dev; - if(estimator->value) MEM_FREE(dev->allocator, estimator->value); - if(estimator->square_value) MEM_FREE(dev->allocator, estimator->value); + if(estimator->value) + MEM_FREE(dev->allocator, estimator->value); + if(estimator->square_value) + MEM_FREE(dev->allocator, estimator->square_value); + if(estimator->status.expected_value) + MEM_FREE(dev->allocator, estimator->status.expected_value); + if(estimator->status.variance) + MEM_FREE(dev->allocator, estimator->status.variance); + if(estimator->status.standard_error) + MEM_FREE(dev->allocator, estimator->status.standard_error); MEM_FREE(dev->allocator, estimator); SMC(device_ref_put(dev)); } @@ -80,7 +91,7 @@ smc_solve void* val = NULL; res_T res = RES_OK; - if(!dev || !integrand || !type || !ctx || !out_estimator || !check_type(type)) { + if(!dev || !integrand || !type || !out_estimator || !check_type(type)) { res = RES_BAD_ARG; goto error; } @@ -93,25 +104,27 @@ smc_solve SMC(device_ref_get(dev)); estimator->dev = dev; ref_init(&estimator->ref); - estimator->value = type->create(dev->allocator); - if(!estimator->value) { - res = RES_MEM_ERR; - goto error; - } - estimator->square_value = type->create(dev->allocator); - if(!estimator->value) { - res = RES_MEM_ERR; - goto error; - } - estimator->nsamples = 0; - val = type->create(dev->allocator); - if(!val) { - res = RES_MEM_ERR; - goto error; - } + #define TYPE_CREATE(Dst) { \ + (Dst) = type->create(dev->allocator); \ + if(!(Dst)) { \ + res = RES_MEM_ERR; \ + goto error; \ + } \ + type->zero((Dst)); \ + } (void)0 + TYPE_CREATE(estimator->value); + TYPE_CREATE(estimator->square_value); + TYPE_CREATE(estimator->status.expected_value); + TYPE_CREATE(estimator->status.variance); + TYPE_CREATE(estimator->status.standard_error); + TYPE_CREATE(val); + #undef TYPE_CREATE + estimator->nsamples = 0; + estimator->status.samples_count = 0; + estimator->type = *type; - FOR_EACH(i, 0, 32) { + FOR_EACH(i, 0, 4096) { if(estimator->nsamples == ULONG_MAX) break; integrand(val, ctx); type->add(estimator->value, estimator->value, val); @@ -119,6 +132,7 @@ smc_solve type->add(estimator->square_value, estimator->square_value, val); ++estimator->nsamples; } + type->destroy(dev->allocator, val); exit: if(out_estimator) *out_estimator = estimator; @@ -147,3 +161,46 @@ smc_estimator_ref_put(struct smc_estimator* estimator) return RES_OK; } +res_T +smc_estimator_get_status + (struct smc_estimator* estimator, + struct smc_estimator_status* status) +{ + if(!estimator || !status) + return RES_BAD_ARG; + + if(estimator->nsamples != estimator->status.samples_count) { + estimator->status.samples_count = estimator->nsamples; + /* Variance */ + estimator->type.divi + (estimator->status.expected_value, + estimator->square_value, + estimator->nsamples); + estimator->type.mul + (estimator->status.variance, estimator->value, estimator->value); + estimator->type.divi + (estimator->status.variance, + estimator->status.variance, + estimator->nsamples * estimator->nsamples); + estimator->type.sub + (estimator->status.variance, + estimator->status.expected_value, + estimator->status.variance); + /* Standard error */ + estimator->type.divi + (estimator->status.standard_error, + estimator->status.variance, + estimator->nsamples); + estimator->type.sqrt + (estimator->status.standard_error, + estimator->status.standard_error); + /* Expected value */ + estimator->type.divi + (estimator->status.expected_value, + estimator->value, + estimator->nsamples); + } + *status = estimator->status; + return RES_OK; +} + diff --git a/src/smc_type.c b/src/smc_type.c @@ -31,6 +31,7 @@ #include "smc.h" #include <rsys/mem_allocator.h> +#include <math.h> /******************************************************************************* * SMC float definition @@ -91,6 +92,13 @@ smc_float_divi(void* result, const void* op0, const unsigned long op1) *(float*)result = (float)((double)(*(const float*)op0) / (double)op1); } +static void +smc_float_sqrt(void* result, const void* value) +{ + ASSERT(result && value && *(const float*)value >= 0.f); + *(float*)result = (float)sqrt(*(const float*)value); +} + /******************************************************************************* * Exported constant ******************************************************************************/ @@ -102,6 +110,7 @@ const struct smc_type smc_float = { smc_float_add, smc_float_sub, smc_float_mul, - smc_float_divi + smc_float_divi, + smc_float_sqrt }; diff --git a/src/smc_type_c.h b/src/smc_type_c.h @@ -44,7 +44,8 @@ check_type(const struct smc_type* type) && type->add != NULL && type->sub != NULL && type->mul != NULL - && type->divi != NULL; + && type->divi != NULL + && type->sqrt != NULL; } #endif /* SMC_TYPE_C_H */ diff --git a/src/test_smc_device.c b/src/test_smc_device.c @@ -33,7 +33,6 @@ #include "test_smc_utils.h" #include <rsys/logger.h> -#include <rsys/mem_allocator.h> static void log_stream(const char* msg, void* ctx)