commit f577a0eb98ec3c3f59f7470e948bcbbc7c6c2125
parent 391b88739003e4be5c069784988e94fbafa5edd1
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 7 Sep 2017 17:11:19 +0200
BugFix: use deferred accumulation for MC weights (struct mc_data).
We cannot add several contributions to a MC weight for a single
realisation. We must instead sum all these contributions in a tmp
and contribute to Sum(x) and Sum(sqr_x) once the realisation is
completed.
Because contrib_1^2 + ... + contrib_n^2 != (contrib_1 + ... +
contrib_n)^2.
Diffstat:
5 files changed, 100 insertions(+), 56 deletions(-)
diff --git a/src/ssol.h b/src/ssol.h
@@ -1103,7 +1103,7 @@ ssol_estimator_ref_put
SSOL_API res_T
ssol_estimator_get_mc_global
- (const struct ssol_estimator* estimator,
+ (struct ssol_estimator* estimator,
struct ssol_mc_global* mc_global);
SSOL_API res_T
diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c
@@ -107,15 +107,17 @@ ssol_estimator_ref_put(struct ssol_estimator* estimator)
res_T
ssol_estimator_get_mc_global
- (const struct ssol_estimator* estimator,
+ (struct ssol_estimator* estimator,
struct ssol_mc_global* global)
{
if(!estimator || !global) return RES_BAD_ARG;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &estimator->Name; \
- global->Name.E = data->weight / N; \
- global->Name.V = data->sqr_weight / N - global->Name.E*global->Name.E; \
+ struct mc_data* data = &estimator->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ global->Name.E = weight / N; \
+ global->Name.V = sqr_weight / N - global->Name.E*global->Name.E; \
global->Name.V = global->Name.V > 0 ? global->Name.V : 0; \
global->Name.SE = sqrt(global->Name.V / N); \
} (void)0
@@ -165,9 +167,11 @@ ssol_estimator_get_mc_sampled_x_receiver
mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &mc_rcv1->Name; \
- rcv->Name.E = data->weight / N; \
- rcv->Name.V = data->sqr_weight / N - rcv->Name.E*rcv->Name.E; \
+ struct mc_data* data = &mc_rcv1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ rcv->Name.E = weight / N; \
+ rcv->Name.V = sqr_weight / N - rcv->Name.E*rcv->Name.E; \
rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \
rcv->Name.SE = sqrt(rcv->Name.V / N); \
} (void)0
@@ -236,9 +240,11 @@ ssol_estimator_get_mc_sampled
sampled->nb_samples = mc->nb_samples;
#define SETUP_MC_RESULT(Name, Count) { \
const double N = (double)(Count); \
- const struct mc_data* data = &mc->Name; \
- sampled->Name.E = data->weight / N; \
- sampled->Name.V = data->sqr_weight/N - sampled->Name.E*sampled->Name.E; \
+ struct mc_data* data = &mc->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ sampled->Name.E = weight / N; \
+ sampled->Name.V = sqr_weight/N - sampled->Name.E*sampled->Name.E; \
sampled->Name.V = sampled->Name.V > 0 ? sampled->Name.V : 0; \
sampled->Name.SE = sqrt(sampled->Name.V / N); \
} (void)0
diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h
@@ -29,10 +29,13 @@ struct ssol_instance;
/* Monte carlo data */
struct mc_data {
- double weight;
- double sqr_weight;
+ size_t irealisation;
+ double tmp;
+ /* internal data; use get() */
+ double __weight;
+ double __sqr_weight;
};
-#define MC_DATA_NULL__ { 0, 0 }
+#define MC_DATA_NULL__ { SIZE_MAX, 0, 0, 0 }
static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__;
#define MC_RECEIVER_DATA \
@@ -60,6 +63,56 @@ static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__;
MC_DATA_NULL__
/*******************************************************************************
+ * Deferred MC data accumulators
+ ******************************************************************************/
+static INLINE void
+mc_data_init(struct mc_data* data)
+{
+ ASSERT(data);
+ *data = MC_DATA_NULL;
+}
+
+static INLINE void
+mc_data_flush(struct mc_data* data)
+{
+ ASSERT(data);
+ data->__weight += data->tmp;
+ data->__sqr_weight += data->tmp * data->tmp;
+ data->tmp = 0;
+}
+
+static INLINE void
+mc_data_add_weight(struct mc_data* data, size_t irealisation, double w)
+{
+ ASSERT(data);
+ ASSERT(irealisation != SIZE_MAX);
+ if(irealisation != data->irealisation) {
+ mc_data_flush(data);
+ data->irealisation = irealisation;
+ }
+ data->tmp = w;
+}
+
+static INLINE void
+mc_data_accum(struct mc_data* dst, struct mc_data* src)
+{
+ ASSERT(dst && src);
+ mc_data_flush(dst);
+ mc_data_flush(src);
+ dst->__weight += src->__weight;
+ dst->__sqr_weight += src->__sqr_weight;
+}
+
+static INLINE void
+mc_data_get(struct mc_data* data, double* weight, double* sqr_weight)
+{
+ ASSERT(data && weight && sqr_weight);
+ mc_data_flush(data);
+ *weight = data->__weight;
+ *sqr_weight = data->__sqr_weight;
+}
+
+/*******************************************************************************
* One sided per shape MC data
******************************************************************************/
struct mc_primitive_1side {
diff --git a/src/ssol_mc_receiver.c b/src/ssol_mc_receiver.c
@@ -53,9 +53,11 @@ ssol_estimator_get_mc_receiver
mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &mc_rcv1->Name; \
- rcv->Name.E = data->weight / N; \
- rcv->Name.V = data->sqr_weight/N - rcv->Name.E*rcv->Name.E; \
+ struct mc_data* data = &mc_rcv1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ rcv->Name.E = weight / N; \
+ rcv->Name.V = sqr_weight/N - rcv->Name.E*rcv->Name.E; \
rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \
rcv->Name.SE = sqrt(rcv->Name.V / N); \
} (void)0
@@ -150,9 +152,11 @@ ssol_mc_shape_get_mc_primitive
#define SETUP_MC_RESULT(Name) { \
const double N = (double)shape->N__; \
- const struct mc_data* data = &mc_prim1->Name; \
- prim->Name.E = data->weight / N; \
- prim->Name.V = data->sqr_weight/N - prim->Name.E*prim->Name.E; \
+ struct mc_data* data = &mc_prim1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ prim->Name.E = weight / N; \
+ prim->Name.V = sqr_weight/N - prim->Name.E*prim->Name.E; \
prim->Name.V = prim->Name.V > 0 ? prim->Name.V : 0; \
prim->Name.SE = sqrt(prim->Name.V / N); \
prim->Name.E /= area; \
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -581,10 +581,7 @@ accum_mc_receivers_1side
res_T res = RES_OK;
ASSERT(dst && src);
- #define ACCUM_WEIGHT(Name) { \
- dst->Name.weight += src->Name.weight; \
- dst->Name.sqr_weight += src->Name.sqr_weight; \
- } (void)0
+ #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name)
#define ACCUM_ALL { \
ACCUM_WEIGHT(incoming_flux); \
ACCUM_WEIGHT(incoming_if_no_atm_loss); \
@@ -627,10 +624,8 @@ accum_mc_receivers_1side
res = mc_shape_1side_get_mc_primitive(mc_shape1_dst, iprim, &mc_prim1_dst);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name) { \
- mc_prim1_dst->Name.weight += mc_prim1_src->Name.weight; \
- mc_prim1_dst->Name.sqr_weight += mc_prim1_src->Name.sqr_weight; \
- } (void)0
+ #define ACCUM_WEIGHT(Name)\
+ mc_data_accum(&mc_prim1_dst->Name, &mc_prim1_src->Name)
ACCUM_ALL;
#undef ACCUM_WEIGHT
#undef ACCUM_ALL
@@ -656,10 +651,7 @@ accum_mc_sampled(struct mc_sampled* dst, struct mc_sampled* src)
mc_receiver_init(NULL, &mc_rcv_null);
- #define ACCUM_WEIGHT(Name) { \
- dst->Name.weight += src->Name.weight; \
- dst->Name.sqr_weight += src->Name.sqr_weight; \
- } (void)0
+ #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name)
ACCUM_WEIGHT(cos_factor);
ACCUM_WEIGHT(shadowed);
#undef ACCUM_WEIGHT
@@ -725,10 +717,8 @@ update_mc
res = get_mc_receiver_1side(&thread_ctx->mc_rcvs, pt->inst, pt->side, &mc_rcv1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_rcv1->Name.weight += (W); \
- mc_rcv1->Name.sqr_weight += (W)*(W); \
- } (void)0
+ #define ACCUM_WEIGHT(Name, W)\
+ mc_data_add_weight(&mc_rcv1->Name, irealisation, W)
#define ACCUM_ALL { \
ACCUM_WEIGHT(incoming_flux, pt->incoming_flux); \
ACCUM_WEIGHT(incoming_if_no_atm_loss, pt->incoming_if_no_atm_loss); \
@@ -755,10 +745,8 @@ update_mc
(pt->mc_samp, pt->inst, pt->side, &mc_samp_x_rcv1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_samp_x_rcv1->Name.weight += (W); \
- mc_samp_x_rcv1->Name.sqr_weight += (W)*(W); \
- } (void)0
+ #define ACCUM_WEIGHT(Name, W)\
+ mc_data_add_weight(&mc_samp_x_rcv1->Name, irealisation, W)
ACCUM_ALL;
#undef ACCUM_WEIGHT
@@ -773,10 +761,8 @@ update_mc
res = mc_shape_1side_get_mc_primitive(mc_shape1, pt->prim.prim_id, &mc_prim1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_prim1->Name.weight += (W); \
- mc_prim1->Name.sqr_weight += (W)*(W); \
- } (void)0
+ #define ACCUM_WEIGHT(Name, W)\
+ mc_data_add_weight(&mc_prim1->Name, irealisation, W)
ACCUM_ALL;
#undef ACCUM_WEIGHT
#undef ACCUM_ALL
@@ -790,7 +776,7 @@ error:
static res_T
trace_radiative_path
- (const size_t path_id, /* Unique id of the radiative path */
+ (const size_t irealisation, /* Unique id of the realisation */
const double sampled_area_proxy, /* Overall area of the sampled geometries */
struct thread_context* thread_ctx,
struct ssol_scene* scn,
@@ -837,10 +823,7 @@ trace_radiative_path
if(res != RES_OK) goto error;
}
- #define ACCUM_WEIGHT(Res, W) { \
- Res.weight += (W); \
- Res.sqr_weight += (W)*(W); \
- } (void)0
+ #define ACCUM_WEIGHT(Res, W) mc_data_add_weight(&Res, irealisation, W)
if(!is_lit) { /* The starting point is not lit */
ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.initial_flux);
@@ -889,7 +872,7 @@ trace_radiative_path
/* If receiver update MC results */
if(hit_receiver) {
hit_a_receiver = 1;
- res = update_mc(&pt, path_id, depth, thread_ctx, output);
+ res = update_mc(&pt, irealisation, depth, thread_ctx, output);
if(res != RES_OK) goto error;
} else {
pt.partial_other += pt.incoming_flux * pt.kabs_at_pt;
@@ -1128,12 +1111,10 @@ ssol_solve
/* Merge per thread global MC estimations */
FOR_EACH(i, 0, nthreads) {
- const struct thread_context* thread_ctx;
- thread_ctx = darray_thread_ctx_cdata_get(&thread_ctxs)+i;
- #define ACCUM_WEIGHT(Name) { \
- estimator->Name.weight += thread_ctx->Name.weight; \
- estimator->Name.sqr_weight += thread_ctx->Name.sqr_weight; \
- } (void)0
+ struct thread_context* thread_ctx;
+ thread_ctx = darray_thread_ctx_data_get(&thread_ctxs)+i;
+ #define ACCUM_WEIGHT(Name) \
+ mc_data_accum(&estimator->Name, &thread_ctx->Name)
ACCUM_WEIGHT(cos_factor);
ACCUM_WEIGHT(absorbed_by_receivers);
ACCUM_WEIGHT(shadowed);
