commit 4b1776d27045791c8ebb50170c6f5a7ba02ad4cb
parent 2060dfe39b71247612b73b589649ef0281e95dd9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 3 Jun 2019 16:18:59 +0200
Test the RNG state saved by the estimator
Diffstat:
4 files changed, 67 insertions(+), 5 deletions(-)
diff --git a/src/ssol.h b/src/ssol.h
@@ -1204,7 +1204,7 @@ ssol_mc_shape_get_mc_primitive
SSOL_API res_T
ssol_solve
(struct ssol_scene* scn,
- struct ssp_rng* rng,
+ const struct ssp_rng* rng,
const size_t realisations_count,
const size_t max_failed_count,
const struct ssol_path_tracker* tracker, /* NULL<=>Do not record the paths */
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -1119,7 +1119,7 @@ error:
res_T
ssol_solve
(struct ssol_scene* scn,
- struct ssp_rng* rng_state,
+ const struct ssp_rng* rng_state,
const size_t realisations_count,
const size_t max_failed_count,
const struct ssol_path_tracker* path_tracker,
diff --git a/src/test_ssol_solver1.c b/src/test_ssol_solver1.c
@@ -49,8 +49,13 @@ main(int argc, char** argv)
{
struct spectrum_desc desc = {0};
struct mem_allocator allocator;
+ FILE* stream;
struct ssol_device* dev;
struct ssp_rng* rng;
+ struct ssp_rng* rng2;
+ const struct ssp_rng* rng_state;
+ struct ssp_rng_type rng_type0;
+ struct ssp_rng_type rng_type1;
struct ssol_scene* scene;
struct ssol_shape* dummy;
struct ssol_shape* square;
@@ -237,6 +242,23 @@ main(int argc, char** argv)
CHK(ssol_estimator_get_mc_global(estimator, NULL) == RES_BAD_ARG);
CHK(ssol_estimator_get_mc_global(estimator, &mc_global) == RES_OK);
+ CHK(ssol_estimator_get_rng_state(NULL, &rng_state) == RES_BAD_ARG);
+ CHK(ssol_estimator_get_rng_state(estimator, NULL) == RES_BAD_ARG);
+ CHK(ssol_estimator_get_rng_state(estimator, &rng_state) == RES_OK);
+ CHK(ssp_rng_get_type(rng_state, &rng_type0) == RES_OK);
+ CHK(ssp_rng_get_type(rng, &rng_type1) == RES_OK);
+ CHK(ssp_rng_type_eq(&rng_type0, &rng_type1));
+
+ /* Clone the rng_state */
+ CHK(stream = tmpfile());
+ CHK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng2) == RES_OK);
+ CHK(ssp_rng_write(rng_state, stream) == RES_OK);
+ rewind(stream);
+ CHK(ssp_rng_read(rng2, stream) == RES_OK);
+ CHK(fclose(stream) == 0);
+ CHK(ssp_rng_get(rng2) != ssp_rng_get(rng));
+ CHK(ssp_rng_ref_put(rng2) == RES_OK);
+
CHK(ssol_estimator_ref_get(NULL) == RES_BAD_ARG);
CHK(ssol_estimator_ref_get(estimator) == RES_OK);
CHK(ssol_estimator_ref_put(NULL) == RES_BAD_ARG);
@@ -324,8 +346,8 @@ main(int argc, char** argv)
CHK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv) == RES_OK);
printf("Ir(target) = %g +/- %g\n",
mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
- CHK(eq_eps(mc_rcv.incoming_flux.E, m, 2 * std) == 1);
- CHK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-1) == 1);
+ CHK(eq_eps(mc_rcv.incoming_flux.E, m, 3 * std) == 1);
+ CHK(eq_eps(mc_rcv.incoming_flux.SE, std, std*1e-2) == 1);
CHK(ssol_estimator_ref_put(estimator) == RES_OK);
/* Sample primary mirror only; variance is low */
diff --git a/src/test_ssol_solver6.c b/src/test_ssol_solver6.c
@@ -51,6 +51,7 @@ main(int argc, char** argv)
struct mem_allocator allocator;
struct ssol_device* dev;
struct ssp_rng* rng;
+ const struct ssp_rng* rng_state;
struct ssol_object *m_object1;
struct ssol_object *m_object2;
struct ssol_object* t_object1;
@@ -74,11 +75,14 @@ main(int argc, char** argv)
struct ssol_sun* sun;
struct ssol_spectrum* spectrum;
struct ssol_estimator* estimator;
+ struct ssol_estimator* estimator2;
struct ssol_mc_global mc_global;
+ struct ssol_mc_global mc_global2;
struct ssol_mc_receiver mc_rcv;
double dir[3];
double transform[12]; /* 3x4 column major matrix */
-
+ double sum_w, sum_w2, E, V, SE;
+ size_t count;
(void) argc, (void) argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -177,6 +181,8 @@ main(int argc, char** argv)
#define N__ 10000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
CHK(ssol_solve(scene, rng, N__, 0, NULL, &estimator) == RES_OK);
+ CHK(ssol_estimator_get_realisation_count(estimator, &count) == RES_OK);
+ CHK(count == N__);
CHK(ssol_estimator_get_mc_global(estimator, &mc_global) == RES_OK);
PRINT_GLOBAL(mc_global);
CHK(eq_eps(mc_global.shadowed.E, 100000, 2 * 100000/sqrt(N__)) == 1);
@@ -194,6 +200,39 @@ main(int argc, char** argv)
CHK(eq_eps(mc_rcv.incoming_flux.E, 0, 1) == 1);
CHK(mc_rcv.incoming_flux.E == mc_rcv.absorbed_flux.E);
+ /* Launch another run that is statistically independent of the first one and
+ * that uses 3 more times samples */
+ CHK(ssol_estimator_get_rng_state(estimator, &rng_state) == RES_OK);
+ CHK(ssol_solve(scene, rng_state, N__*3, 0, NULL, &estimator2) == RES_OK);
+
+ /* Check the estimator result */
+ CHK(ssol_estimator_get_realisation_count(estimator2, &count) == RES_OK);
+ CHK(count == N__*3);
+ CHK(ssol_estimator_get_mc_global(estimator2, &mc_global2) == RES_OK);
+ CHK(eq_eps(mc_global2.shadowed.E, 100000, 2 * 100000/sqrt(N__)) == 1);
+ CHK(eq_eps(mc_global.shadowed.SE/sqrt(3), mc_global2.shadowed.SE, 1.e-1));
+
+ CHK(mc_global.shadowed.E != mc_global2.shadowed.E);
+ CHK(mc_global.shadowed.SE != mc_global2.shadowed.SE);
+
+ /* Merge the 2 estimations */
+ sum_w = mc_global.shadowed.E * N__;
+ sum_w += mc_global2.shadowed.E * N__*3;
+ V = mc_global.shadowed.SE * mc_global.shadowed.SE * N__;
+ sum_w2 = (V + mc_global.shadowed.E * mc_global.shadowed.E) * N__;
+ V = mc_global2.shadowed.SE * mc_global2.shadowed.SE * N__*3;
+ sum_w2 += (V + mc_global2.shadowed.E * mc_global2.shadowed.E) * N__*3;
+ E = sum_w / (4*N__);
+ V = sum_w2 / (4*N__) - E*E;
+ SE = sqrt(V/(4*N__));
+
+ /* Check that the 2 runs are effectively independent, i.e. check the
+ * convergence ratio. By combining the 2 runs we have 4 more times samples,
+ * the standard deviation must be thus devided by 2 while the estimate must
+ * be compatible with the right value regarding the new error */
+ CHK(eq_eps(E, 100000, 3*SE));
+ CHK(eq_eps(SE, mc_global.shadowed.SE / 2, SE*1.e-2));
+
/* Free data */
CHK(ssol_instance_ref_put(heliostat1) == RES_OK);
CHK(ssol_instance_ref_put(target1) == RES_OK);
@@ -210,6 +249,7 @@ main(int argc, char** argv)
CHK(ssol_material_ref_put(v_mtl) == RES_OK);
CHK(ssol_device_ref_put(dev) == RES_OK);
CHK(ssol_estimator_ref_put(estimator) == RES_OK);
+ CHK(ssol_estimator_ref_put(estimator2) == RES_OK);
CHK(ssol_scene_ref_put(scene) == RES_OK);
CHK(ssp_rng_ref_put(rng) == RES_OK);
CHK(ssol_spectrum_ref_put(spectrum) == RES_OK);
