stardis-solver

commit fa745a4d975304f19fdebac53383fc99fb96226a
parent 5795f52a0e2e533cf58d0b5720dd24f81fb01547
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 27 Feb 2019 10:42:05 +0100

Fix a possible correlation in sdis_green_function_solve

Diffstat:
M
src/sdis_green.c
 | 
44
+++++++++++++++++++++++++++++++++++++-------
M
src/sdis_green.h
 | 
7
+++++--
M
src/sdis_heat_path_radiative_Xd.h
 | 
2
+-
M
src/sdis_solve_Xd.h
 | 
2
+-
M
src/test_sdis_solve_probe.c
 | 
2
--

5 files changed, 44 insertions(+), 13 deletions(-)
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -191,6 +191,9 @@ struct sdis_green_function {
   size_t npaths_valid;
   size_t npaths_invalid;
 
+  struct ssp_rng_type rng_type;
+  FILE* rng_state;
+
   ref_T ref;
   struct sdis_device* dev;
 };
@@ -408,6 +411,7 @@ green_function_release(ref_T* ref)
   htable_medium_release(&green->media);
   htable_interf_release(&green->interfaces);
   darray_green_path_release(&green->paths);
+  if(green->rng_state) fclose(green->rng_state);
   MEM_RM(dev->allocator, green);
   SDIS(device_ref_put(dev));
 }
@@ -452,9 +456,13 @@ sdis_green_function_solve
     goto error;
   }
 
-  /* FIXME do not use a new RNG. Save the RNG state into the green function
-   * after its estimation and initialize the following rng with this state */
-  res = ssp_rng_create(green->dev->allocator, &ssp_rng_mt19937_64, &rng);
+  res = ssp_rng_create(green->dev->allocator, &green->rng_type, &rng);
+  if(res != RES_OK) goto error;
+
+  /* Avoid correlation by defining the RNG state from the final state of the
+   * RNG used to estimate the green function */
+  rewind(green->rng_state);
+  res = ssp_rng_read(rng, green->rng_state);
   if(res != RES_OK) goto error;
 
   npaths = darray_green_path_size_get(&green->paths);
@@ -464,7 +472,7 @@ sdis_green_function_solve
   if(res != RES_OK) goto error;
 
   /* Solve the green function */
-  FOR_EACH(ipath, 0, npaths) { /* TODO add multi-threading (?) */
+  FOR_EACH(ipath, 0, npaths) {
     const double time = sample_time(rng, time_range);
     double w;
 
@@ -685,6 +693,12 @@ green_function_create
   green->npaths_valid = SIZE_MAX;
   green->npaths_invalid = SIZE_MAX;
 
+  green->rng_state = tmpfile();
+  if(!green->rng_state) {
+    res = RES_IO_ERR;
+    goto error;
+  }
+
 exit:
   *out_green = green;
   return res;
@@ -759,12 +773,24 @@ error:
 }
 
 res_T
-green_function_finalize(struct sdis_green_function* green)
+green_function_finalize
+  (struct sdis_green_function* green,
+   struct ssp_rng_proxy* proxy)
 {
   size_t i, n;
+  res_T res = RES_OK;
 
-  if(!green) return RES_BAD_ARG;
+  if(!green || !proxy) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
 
+  /* Save the RNG state */
+  SSP(rng_proxy_get_type(proxy, &green->rng_type));
+  res = ssp_rng_proxy_write(proxy, green->rng_state);
+  if(res != RES_OK) goto error;
+
+  /* Compute the number of valid/invalid green paths */
   green->npaths_valid = 0;
   n = darray_green_path_size_get(&green->paths);
   FOR_EACH(i, 0, n) {
@@ -772,7 +798,11 @@ green_function_finalize(struct sdis_green_function* green)
     green->npaths_valid += path->limit_type != SDIS_POINT_NONE;
   }
   green->npaths_invalid = n - green->npaths_valid;
-  return RES_OK;
+
+exit:
+  return res;
+error:
+  goto exit;
 }
 
 res_T
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -20,6 +20,7 @@
 
 /* Forward declaration */
 struct sdis_green_function;
+struct ssp_rng_proxy;
 struct green_path;
 
 struct green_path_handle {
@@ -41,10 +42,12 @@ green_function_merge_and_clear
   (struct sdis_green_function* dst,
    struct sdis_green_function* src);
 
-/* Finalize the green function state (e.g.: computes the #paths & #failures) */
+/* Finalize the green function state (e.g.: computes the #paths & #failures,
+ * save the rng state, etc.) */
 extern LOCAL_SYM res_T
 green_function_finalize
-  (struct sdis_green_function* green);
+  (struct sdis_green_function* green,
+   struct ssp_rng_proxy* rng_proxy); /* Proxy RNG used to estimate the function */
 
 extern LOCAL_SYM res_T
 green_function_create_path
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -39,7 +39,7 @@ XD(trace_radiative_path)
    struct XD(temperature)* T)
 {
   /* The radiative random walk is always performed in 3D. In 2D, the geometry
-   * are assumed to be extruded to the infinty along the Z dimension. */
+   * are assumed to be extruded to the infinity along the Z dimension. */
   float N[3] = {0, 0, 0};
   float dir[3] = {0, 0, 0};
   res_T res = RES_OK;
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -305,7 +305,7 @@ XD(solve_probe)
     }
 
     /* Finalize the estimated green */
-    res = green_function_finalize(green);
+    res = green_function_finalize(green, rng_proxy);
     if(res != RES_OK) goto error;
   }
 
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -448,8 +448,6 @@ main(int argc, char** argv)
   BA(sdis_estimator_for_each_path(estimator, NULL, &dump_ctx));
   OK(sdis_estimator_for_each_path(estimator, process_heat_path, &dump_ctx));
 
-  dump_heat_paths(stderr, estimator);
-
   OK(sdis_estimator_ref_put(estimator));
   OK(sdis_scene_ref_put(scn));
   OK(sdis_device_ref_put(dev));