star-mc

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

commit aa09f7a21f99058469f04fdc3255cf32d2ba1018
parent f249eb11ac779fffe40f3d85d5b181a877f26e1c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 22 Jan 2018 15:48:43 +0100

Rm the raw integration, upd the integrator API

Rename the `max_steps' attribute of the integrator as `max_realisations'.
Add the thread identifier as input of the integrand.

Diffstat:
Mcmake/CMakeLists.txt | 2--
Msrc/smc.h | 56++------------------------------------------------------
Dsrc/smc_accumulator.c | 237-------------------------------------------------------------------------------
Msrc/smc_estimator.c | 17++++++++++-------
Msrc/smc_type_c.h | 10----------
Msrc/test_smc_doubleN.c | 5+++--
Msrc/test_smc_errors.c | 22++++++++++------------
Dsrc/test_smc_integrate.c | 182-------------------------------------------------------------------------------
Msrc/test_smc_light_path.c | 6++++--
Msrc/test_smc_solve.c | 8+++++---
10 files changed, 34 insertions(+), 511 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -64,7 +64,6 @@ set(VERSION_PATCH 1) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SMC_FILES_SRC - smc_accumulator.c smc_device.c smc_doubleN.c smc_estimator.c @@ -125,7 +124,6 @@ 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}) new_test(test_smc_errors) if(NOT TARGET Star3D) diff --git a/src/smc.h b/src/smc.h @@ -75,15 +75,6 @@ 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_accum { - 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 }; @@ -95,37 +86,19 @@ struct smc_estimator_status { size_t NF; /* Failed samples count */ }; -struct smc_accumulator_status { - void* value; /* Accumulated value */ - size_t N; /* Samples count */ -}; - struct smc_integrator { - res_T (*integrand)(void* value, struct ssp_rng* rng, void* ctx); + res_T (*integrand)(void* val, struct ssp_rng*, const int ithread, void* ctx); const struct smc_type* type; - size_t max_steps; /* Maximum # of experiments */ + size_t max_realisations; /* Maximum # of realisations */ size_t max_failures; /* Cancel integration past # failed samples */ }; static const struct smc_integrator SMC_INTEGRATOR_NULL = { NULL, NULL, SIZE_MAX, SIZE_MAX }; -#define SMC_AUTO_CANCEL_OFF SIZE_MAX - -struct smc_integrator_accum { - 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, 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 @@ -208,31 +181,6 @@ smc_estimator_get_status (struct smc_estimator* estimator, struct smc_estimator_status* status); -/******************************************************************************* - * Accumulation API - ******************************************************************************/ -/* Raw integration. - * The realisations must be accumulated by the user defined integrand */ -SMC_API res_T -smc_integrate - (struct smc_device* dev, - struct smc_integrator_accum* 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 @@ -1,237 +0,0 @@ -/* 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 "smc_device_c.h" -#include "smc_type_c.h" - -#include <rsys/mem_allocator.h> - -#include <omp.h> - -struct smc_accumulator { - struct smc_type_accum 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_accum* 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_accum(struct smc_integrator_accum* integrator) -{ - ASSERT(integrator); - return integrator->integrand - && integrator->type - && integrator->max_steps - && check_type_accum(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_accum* integrator, - void* ctx, - struct smc_accumulator** out_accum) -{ - struct smc_accumulator* accum = NULL; - int64_t i; - size_t nthreads = 0; - 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; - goto error; - } - nthreads = device_get_threads_count(dev); - - /* Create and initialize per thread memory variables */ - nsamples = MEM_CALLOC(dev->allocator, nthreads, sizeof(size_t)); - if(!nsamples) { - res = RES_MEM_ERR; - goto error; - } - - accums = MEM_CALLOC(dev->allocator, nthreads, sizeof(void*)); - if(!accums) { - res = RES_MEM_ERR; - goto error; - } - - FOR_EACH(i, 0, (int64_t)nthreads) { - 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 shared(cancel, nfailed) - { - const int ithread = omp_get_thread_num(); - #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; - } - } - } - } - } - - /* 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(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_estimator.c b/src/smc_estimator.c @@ -112,7 +112,7 @@ check_integrator(struct smc_integrator* integrator) ASSERT(integrator); return integrator->integrand && integrator->type - && integrator->max_steps + && integrator->max_realisations && check_type(integrator->type); } @@ -202,9 +202,11 @@ smc_solve { const int ithread = omp_get_thread_num(); #pragma omp for schedule(static) - for(i = 0; i < (int64_t)integrator->max_steps; ++i) { + for(i = 0; i < (int64_t)integrator->max_realisations; ++i) { if(!cancel) { - if(integrator->integrand(vals[ithread], dev->rngs[ithread], ctx) == RES_OK) { + const res_T res_local = integrator->integrand + (vals[ithread], dev->rngs[ithread], ithread, ctx); + if(res_local == RES_OK) { /* call succeded */ integrator->type->add(accums[ithread], accums[ithread], vals[ithread]); integrator->type->mul(vals[ithread], vals[ithread], vals[ithread]); @@ -301,11 +303,12 @@ smc_solve_N #pragma omp for schedule(static, 1) for(i = 0; i < (int64_t)count; ++i) { size_t istep; - FOR_EACH(istep, 0, integrator->max_steps) { + FOR_EACH(istep, 0, integrator->max_realisations) { if(!cancel) { - if(integrator->integrand - (vals[ithread], dev->rngs[ithread], (char*)ctx + (size_t)i*sizeof_ctx) - == RES_OK) { + const res_T res_local = integrator->integrand + (vals[ithread], dev->rngs[ithread], ithread, + (char*)ctx + (size_t)i*sizeof_ctx); + if(res_local == RES_OK) { /* call succeded */ integrator->type->add (estimators[i]->value, estimators[i]->value, vals[ithread]); diff --git a/src/smc_type_c.h b/src/smc_type_c.h @@ -49,15 +49,5 @@ check_type(const struct smc_type* type) && type->sqrt != NULL; } -static INLINE char -check_type_accum(const struct smc_type_accum* type) -{ - ASSERT(type); - return type->create != NULL - && type->destroy != NULL - && type->clear != NULL - && type->add != NULL; -} - #endif /* SMC_TYPE_C_H */ diff --git a/src/test_smc_doubleN.c b/src/test_smc_doubleN.c @@ -45,12 +45,13 @@ enum { struct alpha { double value; }; static res_T -cos_x(void* value, struct ssp_rng* rng, void* context) +cos_x(void* value, struct ssp_rng* rng, const int ithread, void* context) { double* result = SMC_DOUBLEN(value); double samp; struct smc_doubleN_context* ctx = context; const struct alpha* alpha = ctx->integrand_data; + (void)ithread; CHK(value != NULL); CHK(rng != NULL); samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0); @@ -81,7 +82,7 @@ main(int argc, char** argv) integrator.integrand = cos_x; integrator.type = &smc_doubleN; - integrator.max_steps = 100000; + integrator.max_realisations = 100000; alpha.value = a; ctx.count = WEIGHTS_COUNT; ctx.integrand_data = &alpha; diff --git a/src/test_smc_errors.c b/src/test_smc_errors.c @@ -35,17 +35,17 @@ #include <star/ssp.h> static res_T -compute_1(void* value, struct ssp_rng* rng, void* context) +compute_1(void* value, struct ssp_rng* rng, const int ithread, 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; + CHK(ithread >= 0); + if(ok) { + SMC_DOUBLE(value) = 1; /* weight is 1 */ + } else { + SMC_DOUBLE(value) = 0; /* clear weight */ + } /* a realization has been carried out only if ok */ return ok ? RES_OK : RES_BAD_ARG; } @@ -66,9 +66,7 @@ main(int argc, char** argv) integrator.integrand = &compute_1; integrator.type = &smc_double; - integrator.max_steps = 10000; - /* set auto cancel OFF */ - integrator.max_failures = SMC_AUTO_CANCEL_OFF; + integrator.max_realisations = 10000; SMC(solve(dev, &integrator, NULL, &estimator)); SMC(estimator_get_status(estimator, &status)); @@ -83,13 +81,13 @@ main(int argc, char** argv) SMC(estimator_ref_put(estimator)); /* set auto cancel ON */ - integrator.max_failures = integrator.max_steps / 1000; + integrator.max_failures = integrator.max_realisations / 1000; SMC(solve(dev, &integrator, NULL, &estimator)); SMC(estimator_get_status(estimator, &status)); /* it is expected that solve was auto-canceled */ - CHK(integrator.max_steps != status.N); + CHK(integrator.max_realisations != status.N); CHK(status.NF >= integrator.max_failures + 1); /* result must be 1 with std err = 0 */ diff --git a/src/test_smc_integrate.c b/src/test_smc_integrate.c @@ -1,182 +0,0 @@ -/* Copyright (C) |Meso|Star> 2015-2016 (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) -{ - CHK(allocator != NULL); - CHK(ctx == (void*)0xC0DE); - return MEM_CALLOC(allocator, 2, sizeof(double)); -} - -static void -dbl2_destroy(struct mem_allocator* allocator, void* data) -{ - CHK(allocator != NULL); - CHK(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_accum smc_accum_dbl2 = { - dbl2_create, - dbl2_destroy, - dbl2_clear, - dbl2_add -}; - -static res_T -rcp_x(void* value, struct ssp_rng* rng, void* ctx) -{ - double* result = value; - double samp; - double val; - CHK(value != NULL); - CHK(rng != NULL); - CHK(ctx == (void*)0xC0DE); - samp = ssp_rng_uniform_double(rng, 1.0, 4.0); - val = 1.0 / samp * 3; - result[0] += val; - result[1] += val*val; - return RES_OK; -} - -static res_T -cos_x(void* value, struct ssp_rng* rng, void* ctx) -{ - double* result = value; - double samp; - double val; - CHK(value != NULL); - CHK(rng != NULL); - CHK(ctx == (void*)0xC0DE); - samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0); - val = cos(samp) * PI / 2.0; - result[0] += val; - result[1] += val*val; - return RES_OK; -} - -int -main(int argc, char** argv) -{ - struct mem_allocator allocator; - struct smc_device* dev; - struct smc_integrator_accum integrator = SMC_INTEGRATOR_ACCUM_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); - CHK(smc_device_create(NULL, &allocator, SMC_NTHREADS_DEFAULT, NULL, &dev) - == RES_OK); - - integrator.integrand = rcp_x; - integrator.type = &smc_accum_dbl2; - integrator.max_steps = 10000; - - CHK(smc_integrate(NULL, NULL, (void*)0xC0DE, NULL) == RES_BAD_ARG); - CHK(smc_integrate(dev, NULL, (void*)0xC0DE, NULL) == RES_BAD_ARG); - CHK(smc_integrate(NULL, &integrator, (void*)0xC0DE, NULL) == RES_BAD_ARG); - CHK(smc_integrate(dev, &integrator, (void*)0xC0DE, NULL) == RES_BAD_ARG); - CHK(smc_integrate(NULL, NULL, (void*)0xC0DE, &accum) == RES_BAD_ARG); - CHK(smc_integrate(dev, NULL, (void*)0xC0DE, &accum) == RES_BAD_ARG); - CHK(smc_integrate(NULL, &integrator, (void*)0xC0DE, &accum) == RES_BAD_ARG); - CHK(smc_integrate(dev, &integrator, (void*)0xC0DE, &accum) == RES_OK); - - CHK(smc_accumulator_get_status(NULL, NULL) == RES_BAD_ARG); - CHK(smc_accumulator_get_status(accum, NULL) == RES_BAD_ARG); - CHK(smc_accumulator_get_status(NULL, &status) == RES_BAD_ARG); - CHK(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); - CHK(eq_eps(log(4.0) - log(1.0), E, 2.0 * SE) == 1); - - CHK(smc_accumulator_ref_get(NULL) == RES_BAD_ARG); - CHK(smc_accumulator_ref_get(accum) == RES_OK); - CHK(smc_accumulator_ref_put(NULL) == RES_BAD_ARG); - CHK(smc_accumulator_ref_put(accum) == RES_OK); - CHK(smc_accumulator_ref_put(accum) == RES_OK); - - integrator.type = NULL; - CHK(smc_integrate(dev, &integrator, NULL, &accum) == RES_BAD_ARG); - integrator.type = &smc_accum_dbl2; - integrator.integrand = NULL; - CHK(smc_integrate(dev, &integrator, NULL, &accum) == RES_BAD_ARG); - - integrator.integrand = cos_x; - CHK(smc_integrate(dev, &integrator, (void*)0xC0DE, &accum) == RES_OK); - CHK(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); - CHK(eq_eps(2*sin(PI/4.0), E, 2 * SE) == 1); - - CHK(smc_accumulator_ref_put(accum) == RES_OK); - CHK(smc_device_ref_put(dev) == RES_OK); - - check_memory_allocator(&allocator); - mem_shutdown_proxy_allocator(&allocator); - CHK(mem_allocated_size() == 0); - return 0; -} - diff --git a/src/test_smc_light_path.c b/src/test_smc_light_path.c @@ -245,7 +245,8 @@ skydome_lighting(const float wi[3]) } static res_T -light_path_integrator(void* value, struct ssp_rng* rng, void* data) +light_path_integrator + (void* value, struct ssp_rng* rng, const int ithread, void* data) { struct integrator_context* ctx = data; float pix_lower[2], pix_upper[2]; /* Pixel AABB */ @@ -257,6 +258,7 @@ light_path_integrator(void* value, struct ssp_rng* rng, void* data) CHK(value != NULL); CHK(rng != NULL); + CHK(ithread >= 0); CHK(data != NULL); /* Compute the pixel bound */ @@ -375,7 +377,7 @@ main(int argc, char** argv) integrator.integrand = light_path_integrator; integrator.type = &smc_float; - integrator.max_steps = 4; + integrator.max_realisations = 4; contexts = MEM_CALLOC (&allocator, IMG_WIDTH*IMG_HEIGHT, sizeof(struct integrator_context)); diff --git a/src/test_smc_solve.c b/src/test_smc_solve.c @@ -35,12 +35,13 @@ #include <math.h> static res_T -rcp_x(void* value, struct ssp_rng* rng, void* ctx) +rcp_x(void* value, struct ssp_rng* rng, const int ithread, void* ctx) { double* result = value; double samp; CHK(value != NULL); CHK(rng != NULL); + CHK(ithread >= 0); CHK(ctx == NULL); samp = ssp_rng_uniform_double(rng, 1.0, 4.0); *result = 1.0 / samp * 3; @@ -48,12 +49,13 @@ rcp_x(void* value, struct ssp_rng* rng, void* ctx) } static res_T -cos_x(void* value, struct ssp_rng* rng, void* ctx) +cos_x(void* value, struct ssp_rng* rng, const int ithread, void* ctx) { float* result = value; double samp; CHK(value != NULL); CHK(rng != NULL); + CHK(ithread >= 0); CHK(ctx == (void*)0xC0DE); samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0); *result = (float)(cos(samp) * PI / 2.0); @@ -76,7 +78,7 @@ main(int argc, char** argv) integrator.integrand = rcp_x; integrator.type = &smc_double; - integrator.max_steps = 10000; + integrator.max_realisations = 10000; CHK(smc_solve(NULL, NULL, NULL, NULL) == RES_BAD_ARG); CHK(smc_solve(dev, NULL, NULL, NULL) == RES_BAD_ARG);