commit c5ef277973c439a7698525ca44e380eb2f532b38
parent 08be5cd96482f20401945dc96fe0ce807435ca70
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 13 Apr 2016 17:47:16 +0200
API break: the realization function has now to return a status OK/ERROR.
Use it to discard realizations in error.
The failed realizations are not retried.
Estimators now include the number of failed realizations, and integration
has an auto cancel feature based on the number of failures.
Diffstat:
9 files changed, 222 insertions(+), 38 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -126,6 +126,7 @@ if(NOT NO_TEST)
new_test(test_smc_solve ${MATH_LIB})
new_test(test_smc_doubleN ${MATH_LIB})
new_test(test_smc_integrate ${MATH_LIB})
+ new_test(test_smc_errors)
if(NOT TARGET Star3D)
rcmake_copy_runtime_libraries(test_smc_solve)
diff --git a/src/smc.h b/src/smc.h
@@ -91,7 +91,8 @@ struct smc_estimator_status {
void* E; /* Expected value */
void* V; /* Variance */
void* SE; /* Standard error, i.e. sqrt(V / N) */
- size_t N; /* Samples count */
+ size_t N; /* OK samples count */
+ size_t NF; /* Failed samples count */
};
struct smc_accumulator_status {
@@ -100,22 +101,26 @@ struct smc_accumulator_status {
};
struct smc_integrator {
- void (*integrand)(void* value, struct ssp_rng* rng, void* ctx);
+ res_T (*integrand)(void* value, struct ssp_rng* rng, void* ctx);
const struct smc_type* type;
size_t max_steps; /* Maximum # of experiments */
+ size_t max_failures; /* Cancel integration past # failed samples */
};
static const struct smc_integrator SMC_INTEGRATOR_NULL =
-{ NULL, NULL, SIZE_MAX };
+{ NULL, NULL, SIZE_MAX, SIZE_MAX };
+
+#define SMC_AUTO_CANCEL_OFF SIZE_MAX
struct smc_integrator_accum {
- void (*integrand)(void* accum, struct ssp_rng* rng, void* ctx);
+ res_T (*integrand)(void* accum, struct ssp_rng* rng, void* ctx);
const struct smc_type_accum* type;
size_t max_steps; /* Maximum # of accumulation */
+ size_t max_failures; /* Cancel integration past # failed samples */
};
static const struct smc_integrator_accum SMC_INTEGRATOR_ACCUM_NULL =
-{ NULL, NULL, SIZE_MAX };
+{ NULL, NULL, SIZE_MAX, SIZE_MAX };
/* Forward declaration of opaque types */
struct smc_device; /* Entry point of the library */
@@ -174,6 +179,8 @@ smc_device_get_rng_type
/*******************************************************************************
* Integration API
******************************************************************************/
+
+/* Return RES_BAD_OP if auto cancel is enabled and the number of error is reached */
SMC_API res_T
smc_solve
(struct smc_device* dev,
@@ -181,6 +188,7 @@ smc_solve
void* ctx,
struct smc_estimator** estimator);
+/* Return RES_BAD_OP if auto cancel is enabled and the number of error is reached */
SMC_API res_T
smc_solve_N
(struct smc_device* dev,
@@ -206,8 +214,10 @@ smc_estimator_get_status
/*******************************************************************************
* Accumulation API
******************************************************************************/
-/* Raw integration. The realisations must be accumulated by the user defined
- * integrand */
+/* Raw integration.
+ * The realisations must be accumulated by the user defined integrand
+ * Return RES_BAD_OP if auto cancel is enabled and the number of error is reached
+ */
SMC_API res_T
smc_integrate
(struct smc_device* dev,
diff --git a/src/smc_accumulator.c b/src/smc_accumulator.c
@@ -129,6 +129,8 @@ smc_integrate
void** accums = NULL;
size_t* nsamples = NULL;
res_T res = RES_OK;
+ size_t nfailed = 0;
+ int cancel = 0;
if(!dev || !integrator || !out_accum || !check_integrator_accum(integrator)) {
res = RES_BAD_ARG;
@@ -136,7 +138,7 @@ smc_integrate
}
nthreads = device_get_threads_count(dev);
- /* Create and initializdd per thread memory variables */
+ /* Create and initialize per thread memory variables */
nsamples = MEM_CALLOC(dev->allocator, nthreads, sizeof(size_t));
if(!nsamples) {
res = RES_MEM_ERR;
@@ -163,11 +165,26 @@ smc_integrate
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) {
+ #pragma omp parallel shared(cancel, nfailed)
+ {
const int ithread = omp_get_thread_num();
- integrator->integrand(accums[ithread], dev->rngs[ithread], ctx);
- ++nsamples[ithread];
+ #pragma omp for schedule(static)
+ for(i=0; i < (int64_t)integrator->max_steps; ++i) {
+ if(!cancel) {
+ if(integrator->integrand(accums[ithread], dev->rngs[ithread], ctx) == RES_OK) {
+ /* call succeded */
+ ++nsamples[ithread];
+ }
+ else {
+#pragma omp atomic
+ ++nfailed;
+ if (nfailed > integrator->max_failures) {
+ cancel = 1;
+ res = RES_BAD_OP;
+ }
+ }
+ }
+ }
}
/* Merge the parallel estimation into the final estimator */
diff --git a/src/smc_estimator.c b/src/smc_estimator.c
@@ -44,6 +44,7 @@ struct smc_estimator {
void* value;
void* square_value;
size_t nsamples;
+ size_t nfailed;
struct smc_estimator_status status;
@@ -89,7 +90,9 @@ estimator_create
TYPE_CREATE(estimator->status.SE);
#undef TYPE_CREATE
estimator->nsamples = 0;
+ estimator->nfailed = 0;
estimator->status.N = 0;
+ estimator->status.NF = 0;
estimator->type = *type;
exit:
@@ -155,6 +158,8 @@ smc_solve
void** accums_sqr = NULL;
size_t* nsamples = NULL;
res_T res = RES_OK;
+ size_t nfailed = 0;
+ int cancel = 0;
if(!dev || !integrator || !out_estimator || !check_integrator(integrator)) {
res = RES_BAD_ARG;
@@ -193,14 +198,29 @@ smc_solve
if(res != RES_OK) goto error;
/* Parallel evaluation of the simulation */
- #pragma omp parallel for schedule(static)
- for(i = 0; i < (int64_t)integrator->max_steps; ++i) {
+ #pragma omp parallel shared(cancel, nfailed)
+ {
const int ithread = omp_get_thread_num();
- integrator->integrand(vals[ithread], dev->rngs[ithread], ctx);
- integrator->type->add(accums[ithread], accums[ithread], vals[ithread]);
- integrator->type->mul(vals[ithread], vals[ithread], vals[ithread]);
- integrator->type->add(accums_sqr[ithread], accums_sqr[ithread], vals[ithread]);
- ++nsamples[ithread];
+ #pragma omp for schedule(static)
+ for(i = 0; i < (int64_t)integrator->max_steps; ++i) {
+ if (! cancel) {
+ if(integrator->integrand(vals[ithread], dev->rngs[ithread], ctx) == RES_OK) {
+ /* call succeded */
+ integrator->type->add(accums[ithread], accums[ithread], vals[ithread]);
+ integrator->type->mul(vals[ithread], vals[ithread], vals[ithread]);
+ integrator->type->add(accums_sqr[ithread], accums_sqr[ithread], vals[ithread]);
+ ++nsamples[ithread];
+ }
+ else {
+#pragma omp atomic
+ ++nfailed;
+ if (nfailed > integrator->max_failures) {
+ cancel = 1;
+ res = RES_BAD_OP;
+ }
+ }
+ }
+ }
}
/* Merge the parallel estimation into the final estimator */
@@ -210,6 +230,7 @@ smc_solve
integrator->type->add
(estimator->square_value, estimator->square_value, accums_sqr[i]);
}
+ estimator->nfailed = nfailed;
exit:
if(raw) { /* Release temporary variables */
@@ -243,6 +264,7 @@ smc_solve_N
int64_t i;
size_t nthreads = 0;
void** vals = NULL;
+ int cancel = 0;
if(!estimators) {
res = RES_BAD_ARG;
@@ -275,19 +297,35 @@ smc_solve_N
}
/* Parallel estimation of N simulations */
- #pragma omp parallel for schedule(static, 1)
- for(i = 0; i < (int64_t)count; ++i) {
- size_t istep;
+ #pragma omp parallel shared(cancel)
+ {
const int ithread = omp_get_thread_num();
- FOR_EACH(istep, 0, integrator->max_steps) {
- integrator->integrand
- (vals[ithread], dev->rngs[ithread], (char*)ctx + (size_t)i*sizeof_ctx);
- integrator->type->add
- (estimators[i]->value, estimators[i]->value, vals[ithread]);
- integrator->type->mul(vals[ithread], vals[ithread], vals[ithread]);
- integrator->type->add
- (estimators[i]->square_value, estimators[i]->square_value, vals[ithread]);
- ++estimators[i]->nsamples;
+ #pragma omp for schedule(static, 1)
+ for(i = 0; i < (int64_t)count; ++i) {
+ size_t istep;
+ FOR_EACH(istep, 0, integrator->max_steps) {
+ if (! cancel) {
+ if(integrator->integrand
+ (vals[ithread], dev->rngs[ithread], (char*)ctx + (size_t)i*sizeof_ctx)
+ == RES_OK) {
+ /* call succeded */
+ integrator->type->add
+ (estimators[i]->value, estimators[i]->value, vals[ithread]);
+ integrator->type->mul(vals[ithread], vals[ithread], vals[ithread]);
+ integrator->type->add
+ (estimators[i]->square_value, estimators[i]->square_value, vals[ithread]);
+ ++estimators[i]->nsamples;
+ }
+ else {
+#pragma omp atomic
+ ++estimators[i]->nfailed;
+ if (estimators[i]->nfailed > integrator->max_failures) {
+ cancel = 1;
+ res = RES_BAD_OP;
+ }
+ }
+ }
+ }
}
}
@@ -335,8 +373,11 @@ smc_estimator_get_status
if(!estimator || !status)
return RES_BAD_ARG;
- if(estimator->nsamples != estimator->status.N) {
+ if(estimator->nsamples != estimator->status.N
+ || estimator->nfailed != estimator->status.NF)
+ {
estimator->status.N = estimator->nsamples;
+ estimator->status.NF = estimator->nfailed;
/* Variance */
estimator->type.divi
(estimator->status.E, estimator->square_value, estimator->nsamples);
diff --git a/src/test_smc_doubleN.c b/src/test_smc_doubleN.c
@@ -44,7 +44,7 @@ enum {
struct alpha { double value; };
-static void
+static res_T
cos_x(void* value, struct ssp_rng* rng, void* context)
{
double* result = SMC_DOUBLEN(value);
@@ -59,6 +59,8 @@ cos_x(void* value, struct ssp_rng* rng, void* context)
cos(alpha->value*samp) * PI / 2.0;
result[WEIGHT_SENSIBILITY_ALPHA] =
-samp * sin(alpha->value*samp) * PI/2.0;
+
+ return RES_OK;
}
int
diff --git a/src/test_smc_errors.c b/src/test_smc_errors.c
@@ -0,0 +1,108 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (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>
+
+static res_T
+compute_1(void* value, struct ssp_rng* rng, void* context)
+{
+ /* pretend that 10% of the realizations fail */
+ int ok = ssp_rng_canonical(rng) > 0.1;
+ (void)context;
+ if (ok)
+ /* weight is 1 */
+ SMC_DOUBLE(value) = 1;
+ else
+ /* clear weight */
+ SMC_DOUBLE(value) = 0;
+ /* a realization has been carried out only if ok */
+ return ok ? RES_OK : RES_BAD_ARG;
+}
+
+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;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ SMC(device_create(NULL, &allocator, SMC_NTHREADS_DEFAULT, NULL, &dev));
+
+ integrator.integrand = &compute_1;
+ integrator.type = &smc_double;
+ integrator.max_steps = 10000;
+ /* set auto cancel OFF */
+ integrator.max_failures = SMC_AUTO_CANCEL_OFF;
+
+ SMC(solve(dev, &integrator, NULL, &estimator));
+ SMC(estimator_get_status(estimator, &status));
+
+ /* result must be 1 with std err = 0 */
+ CHECK(SMC_DOUBLE(status.E), 1);
+ CHECK(SMC_DOUBLE(status.SE), 0);
+
+ printf("OK realizations = %lu; KO realizations = %lu\n",
+ (unsigned long)status.N, (unsigned long)status.NF);
+
+ SMC(estimator_ref_put(estimator));
+
+ /* set auto cancel ON */
+ integrator.max_failures = integrator.max_steps / 1000;
+
+ CHECK(smc_solve(dev, &integrator, NULL, &estimator),
+ RES_BAD_OP); /* it is expected that solve was canceled */
+ SMC(estimator_get_status(estimator, &status));
+
+ /* result must be 1 with std err = 0 */
+ CHECK(SMC_DOUBLE(status.E), 1);
+ CHECK(SMC_DOUBLE(status.SE), 0);
+ CHECK(status.NF >= integrator.max_failures, 1);
+
+ printf("OK realizations = %lu; KO realizations = %lu\n",
+ (unsigned long)status.N, (unsigned long)status.NF);
+
+ SMC(device_ref_put(dev));
+ SMC(estimator_ref_put(estimator));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
+
diff --git a/src/test_smc_integrate.c b/src/test_smc_integrate.c
@@ -71,7 +71,7 @@ static const struct smc_type_accum smc_accum_dbl2 = {
dbl2_add
};
-static void
+static res_T
rcp_x(void* value, struct ssp_rng* rng, void* ctx)
{
double* result = value;
@@ -84,9 +84,10 @@ rcp_x(void* value, struct ssp_rng* rng, void* ctx)
val = 1.0 / samp * 3;
result[0] += val;
result[1] += val*val;
+ return RES_OK;
}
-static void
+static res_T
cos_x(void* value, struct ssp_rng* rng, void* ctx)
{
double* result = value;
@@ -99,6 +100,7 @@ cos_x(void* value, struct ssp_rng* rng, void* ctx)
val = cos(samp) * PI / 2.0;
result[0] += val;
result[1] += val*val;
+ return RES_OK;
}
int
diff --git a/src/test_smc_light_path.c b/src/test_smc_light_path.c
@@ -243,7 +243,7 @@ skydome_lighting(const float wi[3])
return u * sky_irradiance + (1.f - u) * ground_irradiance;
}
-static void
+static res_T
light_path_integrator(void* value, struct ssp_rng* rng, void* data)
{
struct integrator_context* ctx = data;
@@ -306,6 +306,7 @@ light_path_integrator(void* value, struct ssp_rng* rng, void* data)
throughput *= (float)(ALBEDO / PI) / pdf * cos_theta;
}
*(float*)value = L;
+ return RES_OK;
}
/*******************************************************************************
diff --git a/src/test_smc_solve.c b/src/test_smc_solve.c
@@ -34,7 +34,7 @@
#include <math.h>
-static void
+static res_T
rcp_x(void* value, struct ssp_rng* rng, void* ctx)
{
double* result = value;
@@ -44,9 +44,10 @@ rcp_x(void* value, struct ssp_rng* rng, void* ctx)
CHECK(ctx, NULL);
samp = ssp_rng_uniform_double(rng, 1.0, 4.0);
*result = 1.0 / samp * 3;
+ return RES_OK;
}
-static void
+static res_T
cos_x(void* value, struct ssp_rng* rng, void* ctx)
{
float* result = value;
@@ -56,6 +57,7 @@ cos_x(void* value, struct ssp_rng* rng, void* ctx)
CHECK(ctx, (void*)0xC0DE);
samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0);
*result = (float)(cos(samp) * PI / 2.0);
+ return RES_OK;
}
int
