star-mc

commit f3601b600e7373cdae0d0530b881d3456a15816b
parent afb54ecbd8d3dc1e342b880a753e17c0a2a509ab
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sat, 12 Sep 2015 15:43:11 +0200

Simplify the raw integration process

The accumulation of the realisations is no more automatically performed
by the library, i.e. the user defined integrand must explicitly
accumulate the results in the provided accum value. This may drastically
improve the performances for sparse accum value.

Diffstat:
M
src/smc.h
 | 
2
++
M
src/smc_accumulator.c
 | 
25
+++++++++++--------------
M
src/test_smc_integrate.c
 | 
12
++++++++----

3 files changed, 21 insertions(+), 18 deletions(-)
diff --git a/src/smc.h b/src/smc.h
@@ -191,6 +191,8 @@ smc_estimator_get_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,
diff --git a/src/smc_accumulator.c b/src/smc_accumulator.c
@@ -126,7 +126,6 @@ smc_integrate
   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;
@@ -143,14 +142,19 @@ smc_integrate
     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*));
+  if(!accums) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+
   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; }
+    if(!accums[i]) {
+      res = RES_MEM_ERR;
+      goto error;
+    }
     integrator->type->clear(accums[i]);
   }
 
@@ -162,9 +166,7 @@ smc_integrate
   #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[ithread]);
-    integrator->integrand(vals[ithread], dev->rngs[ithread], ctx);
-    integrator->type->add(accums[ithread], accums[ithread], vals[ithread]);
+    integrator->integrand(accums[ithread], dev->rngs[ithread], ctx);
     ++nsamples[ithread];
   }
 
@@ -176,11 +178,6 @@ smc_integrate
 
 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]);
diff --git a/src/test_smc_integrate.c b/src/test_smc_integrate.c
@@ -76,12 +76,14 @@ rcp_x(void* value, struct ssp_rng* rng, void* ctx)
 {
   double* result = value;
   double samp;
+  double val;
   NCHECK(value, NULL);
   NCHECK(rng, NULL);
   CHECK(ctx, (void*)0xC0DE);
   samp = ssp_rng_uniform_double(rng, 1.0, 4.0);
-  result[0] = 1.0 / samp * 3;
-  result[1] = result[0]*result[0];
+  val =  1.0 / samp * 3;
+  result[0] += val;
+  result[1] += val*val;
 }
 
 static void
@@ -89,12 +91,14 @@ cos_x(void* value, struct ssp_rng* rng, void* ctx)
 {
   double* result = value;
   double samp;
+  double val;
   NCHECK(value, NULL);
   NCHECK(rng, NULL);
   CHECK(ctx, (void*)0xC0DE);
   samp = ssp_rng_uniform_double(rng, -PI/4.0, PI/4.0);
-  result[0] = cos(samp) * PI / 2.0;
-  result[1] = result[0]*result[0];
+  val = cos(samp) * PI / 2.0;
+  result[0] += val;
+  result[1] += val*val;
 }
 
 int