star-mc

commit 615a7b2350a9cbb80d63cac057fd8ef2bc02d611
parent b76078b7f01162598a5901aaf02f80ed1b90697f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  2 Apr 2015 12:05:52 +0200

Provide a Star Sampling RNG to the user defined integrand

Diffstat:
M
cmake/CMakeLists.txt
 | 
5
+++--
M
src/smc.h
 | 
3
++-
M
src/smc_device.c
 | 
6
++++++
M
src/smc_device_c.h
 | 
1
+
M
src/smc_integrator.c
 | 
4
++--
M
src/test_smc_solve.c
 | 
15
+++++++++------

6 files changed, 23 insertions(+), 11 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -41,10 +41,11 @@ option(NO_TEST "Disable the test" OFF)
 ################################################################################
 find_package(RCMake 0.1 REQUIRED)
 find_package(RSys 0.1.1 REQUIRED)
+find_package(StarSP REQUIRED)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
 include(rcmake)
 
-include_directories(${Embree_INCLUDE_DIR} ${RSys_INCLUDE_DIR})
+include_directories(${RSys_INCLUDE_DIR} ${StarSP_INCLUDE_DIR})
 
 ################################################################################
 # Configure and define targets
@@ -75,7 +76,7 @@ set_target_properties(smc PROPERTIES
   VERSION ${VERSION}
   SOVERSION ${VERSION_MAJOR})
 
-target_link_libraries(smc RSys)
+target_link_libraries(smc RSys StarSP)
 if(CMAKE_COMPILER_IS_GNUCC)
   target_link_libraries(smc m)
 endif()
diff --git a/src/smc.h b/src/smc.h
@@ -55,6 +55,7 @@
 /* Forward declaration of external types */
 struct logger;
 struct mem_allocator;
+struct ssp_rng;
 
 /* Generic type descriptor */
 struct smc_type {
@@ -114,7 +115,7 @@ smc_device_ref_put
 SMC_API res_T
 smc_solve
   (struct smc_device* dev,
-   void (*integrand)(void* value, void* ctx),
+   void (*integrand)(void* value, struct ssp_rng* rng, void* ctx),
    const struct smc_type* type,
    void* ctx,
    struct smc_estimator** estimator);
diff --git a/src/smc_device.c b/src/smc_device.c
@@ -32,6 +32,8 @@
 #include "smc.h"
 #include "smc_device_c.h"
 
+#include <star/ssp.h>
+
 #include <rsys/logger.h>
 #include <rsys/mem_allocator.h>
 
@@ -44,6 +46,7 @@ device_release(ref_T* ref)
   struct smc_device* dev;
   ASSERT(ref);
   dev = CONTAINER_OF(ref, struct smc_device, ref);
+  if(dev->rng) ssp_rng_ref_put(dev->rng);
   MEM_FREE(dev->allocator, dev);
 }
 
@@ -73,6 +76,9 @@ smc_device_create
   ref_init(&dev->ref);
   dev->allocator = allocator;
   dev->logger = logger ? logger : LOGGER_DEFAULT;
+  res = ssp_rng_create(allocator, &ssp_rng_mt19937_64, &dev->rng);
+  if(res != RES_OK)
+    goto error;
 
 exit:
   if(out_dev) *out_dev = dev;
diff --git a/src/smc_device_c.h b/src/smc_device_c.h
@@ -35,6 +35,7 @@
 #include <rsys/ref_count.h>
 
 struct smc_device {
+  struct ssp_rng* rng;
   struct logger* logger;
   struct mem_allocator* allocator;
   ref_T ref;
diff --git a/src/smc_integrator.c b/src/smc_integrator.c
@@ -81,7 +81,7 @@ estimator_release(ref_T* ref)
 res_T
 smc_solve
   (struct smc_device* dev,
-   void (*integrand)(void* value, void* ctx),
+   void (*integrand)(void* value, struct ssp_rng* rng, void* ctx),
    const struct smc_type* type,
    void* ctx,
    struct smc_estimator** out_estimator)
@@ -126,7 +126,7 @@ smc_solve
 
   FOR_EACH(i, 0, 4096) {
     if(estimator->nsamples == ULONG_MAX) break;
-    integrand(val, ctx);
+    integrand(val, dev->rng, ctx);
     type->add(estimator->value, estimator->value, val);
     type->mul(val, val, val);
     type->add(estimator->square_value, estimator->square_value, val);
diff --git a/src/test_smc_solve.c b/src/test_smc_solve.c
@@ -32,29 +32,32 @@
 #include "smc.h"
 #include "test_smc_utils.h"
 
+#include <star/ssp.h>
 #include <rsys/math.h>
 
 #include <math.h>
 
 static void
-rcp_x(void* value, void* ctx)
+rcp_x(void* value, struct ssp_rng* rng, void* ctx)
 {
   double* result = value;
-  /* samp in [1, 2] */
-  double samp = ((double)rand() / (double)RAND_MAX) + 1.0;
+  double samp;
+  NCHECK(value, NULL);
+  NCHECK(rng, NULL);
   CHECK(ctx, NULL);
-  (void)ctx;
+  samp = ssp_rng_get_canonical(rng) + 1; /* samp in [1, 2] */
   *result = 1.f / samp;
 }
 
 static void
-cos_x(void* value, void* ctx)
+cos_x(void* value, struct ssp_rng* rng, void* ctx)
 {
   float* result = value;
   /* samp in [PI/4, 3PI/4] */
   float samp = ((float)rand() / (float)RAND_MAX) * (float)(2.0*PI/4.0 + PI/4.0);
+  NCHECK(value, NULL);
+  NCHECK(rng, NULL);
   CHECK(ctx, (void*)0xC0DE);
-  (void)ctx;
   *result = (float)cos(samp);
 }