commit 1dfc0491fc853b9f3c8b9a4e706c3a81bb304eb2
parent 081aa5390949c643ec93587fc9514ee98238c862
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 19 Sep 2017 11:39:25 +0200
BugFix: manage path with unusual depth without bias.
We use a russian roulette to manage both inifine paths with an accurate
0 weight and extra-long paths with an accurate possibly non-0 weight.
The code doesn't rely on some decrease of MC weights over segments, thus
working well with non-physical perfect mirrors.
Diffstat:
2 files changed, 138 insertions(+), 83 deletions(-)
diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h
@@ -71,7 +71,10 @@ mc_data_flush(struct mc_data* data)
}
static INLINE void
-mc_data_add_weight(struct mc_data* data, size_t irealisation, double w)
+mc_data_add_weight
+ (struct mc_data* data,
+ const size_t irealisation,
+ const double w)
{
ASSERT(data);
ASSERT(irealisation != SIZE_MAX);
@@ -102,13 +105,24 @@ mc_data_get(struct mc_data* data, double* weight, double* sqr_weight)
}
static INLINE void
-mc_data_cancel(struct mc_data* data, size_t irealisation)
+mc_data_cancel(struct mc_data* data, const size_t irealisation)
{
ASSERT(data);
if(data->irealisation!=irealisation) return;
data->tmp = 0;
}
+static INLINE void
+mc_data_apply_factor
+ (struct mc_data* data,
+ const size_t irealisation,
+ const double factor)
+{
+ ASSERT(data);
+ if(data->irealisation != irealisation) return;
+ data->tmp *= factor;
+}
+
/*******************************************************************************
* One sided per shape MC data
******************************************************************************/
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -173,6 +173,7 @@ struct point {
const struct ssol_material* material;
/* tmp quantities to compute weights */
double kabs_at_pt;
+ size_t survivor_score;
/* for conservation of energy check */
double energy_loss;
/* MC weights */
@@ -205,7 +206,7 @@ struct point {
{0, 0}, /* UV */ \
0, /* Wavelength */ \
NULL, /* Material */ \
- 0, /* tmp values */ \
+ 0, 0, /* tmp values */ \
0, /* Energy loss */ \
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \
SSOL_FRONT /* Side */ \
@@ -313,6 +314,7 @@ point_init
pt->prev_outgoing_flux = w0;
pt->prev_outgoing_if_no_atm_loss = w0;
pt->prev_outgoing_if_no_field_loss = w0;
+ pt->survivor_score = 0;
d3_set(pt->N, N);
ASSERT(d3_dot(pt->N, pt->dir) <= 0);
@@ -746,6 +748,92 @@ error:
goto exit;
}
+static void
+apply_factor_mc_receiver_1side
+ (struct mc_receiver_1side* rcv,
+ const size_t irealisation,
+ const double factor)
+{
+ mc_data_apply_factor(&rcv->incoming_flux, irealisation, factor);
+ mc_data_apply_factor(&rcv->incoming_if_no_atm_loss, irealisation, factor);
+ mc_data_apply_factor(&rcv->incoming_if_no_field_loss, irealisation, factor);
+ mc_data_apply_factor(&rcv->incoming_lost_in_field, irealisation, factor);
+ mc_data_apply_factor(&rcv->incoming_lost_in_atmosphere, irealisation, factor);
+ mc_data_apply_factor(&rcv->absorbed_flux, irealisation, factor);
+ mc_data_apply_factor(&rcv->absorbed_if_no_atm_loss, irealisation, factor);
+ mc_data_apply_factor(&rcv->absorbed_if_no_field_loss, irealisation, factor);
+ mc_data_apply_factor(&rcv->absorbed_lost_in_field, irealisation, factor);
+ mc_data_apply_factor(&rcv->absorbed_lost_in_atmosphere, irealisation, factor);
+}
+
+static void
+apply_factor_mc
+ (struct thread_context* thread_ctx,
+ const size_t irealisation,
+ const double factor)
+{
+ struct htable_receiver_iterator r_it, r_end;
+ struct htable_sampled_iterator s_it, s_end;
+
+ /* Cancel global MC estimations */
+ mc_data_apply_factor(&thread_ctx->cos_factor, irealisation, factor);
+ mc_data_apply_factor(&thread_ctx->absorbed_by_atmosphere, irealisation, factor);
+ mc_data_apply_factor(&thread_ctx->absorbed_by_receivers, irealisation, factor);
+ mc_data_apply_factor(&thread_ctx->other_absorbed, irealisation, factor);
+ mc_data_apply_factor(&thread_ctx->missing, irealisation, factor);
+ mc_data_apply_factor(&thread_ctx->shadowed, irealisation, factor);
+
+ /* Cancel receiver MC estimations */
+ htable_receiver_begin(&thread_ctx->mc_rcvs, &r_it);
+ htable_receiver_end(&thread_ctx->mc_rcvs, &r_end);
+ while(!htable_receiver_iterator_eq(&r_it, &r_end)) {
+ struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
+ const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
+
+ htable_receiver_iterator_next(&r_it);
+
+ if(inst->receiver_mask & (int)SSOL_FRONT) {
+ apply_factor_mc_receiver_1side(&mc_rcv->front, irealisation, factor);
+ }
+ if(inst->receiver_mask & (int)SSOL_BACK) {
+ apply_factor_mc_receiver_1side(&mc_rcv->back, irealisation, factor);
+ }
+ }
+ /* Cancel sampled instance MC estimations */
+ htable_sampled_begin(&thread_ctx->mc_samps, &s_it);
+ htable_sampled_end(&thread_ctx->mc_samps, &s_end);
+ while(!htable_sampled_iterator_eq(&s_it, &s_end)) {
+ struct mc_sampled* mc_samp = htable_sampled_iterator_data_get(&s_it);
+ htable_sampled_iterator_next(&s_it);
+
+ mc_data_apply_factor(&mc_samp->cos_factor, irealisation, factor);
+ mc_data_apply_factor(&mc_samp->shadowed, irealisation, factor);
+
+ /* dst->by_receiver += src->by_receiver; */
+ htable_receiver_begin(&mc_samp->mc_rcvs, &r_it);
+ htable_receiver_end(&mc_samp->mc_rcvs, &r_end);
+ while(!htable_receiver_iterator_eq(&r_it, &r_end)) {
+ struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
+ const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
+ htable_receiver_iterator_next(&r_it);
+
+ if(inst->receiver_mask & (int)SSOL_FRONT) {
+ apply_factor_mc_receiver_1side(&mc_rcv->front, irealisation, factor);
+ }
+ if(inst->receiver_mask & (int)SSOL_BACK) {
+ apply_factor_mc_receiver_1side(&mc_rcv->back, irealisation, factor);
+ }
+ }
+ }
+}
+static void
+cancel_mc
+ (struct thread_context* thread_ctx,
+ const size_t irealisation)
+{
+ apply_factor_mc(thread_ctx, irealisation, 0);
+}
+
static res_T
trace_radiative_path
(const size_t irealisation, /* Unique id of the realisation */
@@ -765,12 +853,17 @@ trace_radiative_path
struct point pt;
float org[3], dir[3], range[2] = { 0, FLT_MAX };
size_t depth = 0;
+ size_t roulette_interval, typical_max_depth;
int is_lit = 0;
int hit_a_receiver = 0;
+ int killed_by_roulette = 0;
res_T res = RES_OK;
ASSERT(thread_ctx && scn && view_samp && view_rt && ran_sun_dir && ran_sun_wl);
if(tracker) path_init(scn->dev->allocator, &path);
+
+ typical_max_depth = 16; /* This one could come through scn */
+ roulette_interval = 4 * typical_max_depth; /* First roulette */
/* Find a new starting point of the radiative random walk */
res = point_init(&pt, sampled_area_proxy, scn, &thread_ctx->mc_samps,
@@ -928,8 +1021,27 @@ trace_radiative_path
}
depth += !hit_virtual;
- /* FIXME: create a true cancel path for this MC sample */
- ASSERT(depth < 100);
+ if(depth > roulette_interval && depth % roulette_interval == 1) {
+ /* This could be in an infinite reflection path.
+ * To avoid to crash the app while preserving MC weights
+ * we have to use a russian roulette:
+ * 1/2 probability to end the path now compensated by a 2x factor
+ * on weights if the path eventually ends somewhere.
+ * We could have written a more traditional russian roulette that
+ * relies on not applying kabs VS setting weights=0, but this
+ * doesn't work with kabs=0.
+ * The present code works even if weigths remain unchanged
+ * along the path. */
+ double p = ssp_rng_canonical(thread_ctx->rng);
+ if(p > 0.5) {
+ pt.survivor_score += 1; /* This path survived once more */
+ } else {
+ cancel_mc(thread_ctx, irealisation);
+ killed_by_roulette = 1;
+ roulette_interval = typical_max_depth;
+ goto exit; /* break is not enough */
+ }
+ }
/* Update the point */
point_update_from_hit(&pt, scn, org, dir, &hit, &ray_data);
@@ -958,8 +1070,14 @@ trace_radiative_path
/* Check conservation of energy at the realisation level */
ASSERT((double)depth*DBL_EPSILON*pt.initial_flux >= fabs(pt.energy_loss));
+ /* this realisation count account for many that where canceled */
+ if(pt.survivor_score) {
+ const double factor = (double)(1 << pt.survivor_score);
+ apply_factor_mc(thread_ctx, irealisation, factor);
+ }
+
exit:
- if(tracker) {
+ if(tracker && !killed_by_roulette) {
res_T tmp_res = path_register_and_clear(&thread_ctx->paths, &path);
if(tmp_res != RES_OK && res == RES_OK) {
res = tmp_res;
@@ -977,83 +1095,6 @@ error:
goto exit;
}
-static void
-cancel_mc_receiver_1side
- (struct mc_receiver_1side* rcv,
- size_t irealisation)
-{
- mc_data_cancel(&rcv->incoming_flux, irealisation);
- mc_data_cancel(&rcv->incoming_if_no_atm_loss, irealisation);
- mc_data_cancel(&rcv->incoming_if_no_field_loss, irealisation);
- mc_data_cancel(&rcv->incoming_lost_in_field, irealisation);
- mc_data_cancel(&rcv->incoming_lost_in_atmosphere, irealisation);
- mc_data_cancel(&rcv->absorbed_flux, irealisation);
- mc_data_cancel(&rcv->absorbed_if_no_atm_loss, irealisation);
- mc_data_cancel(&rcv->absorbed_if_no_field_loss, irealisation);
- mc_data_cancel(&rcv->absorbed_lost_in_field, irealisation);
- mc_data_cancel(&rcv->absorbed_lost_in_atmosphere, irealisation);
-}
-
-static void
-cancel_mc
- (struct thread_context* thread_ctx,
- size_t irealisation)
-{
- struct htable_receiver_iterator r_it, r_end;
- struct htable_sampled_iterator s_it, s_end;
-
- /* Cancel global MC estimations */
- mc_data_cancel(&thread_ctx->cos_factor, irealisation);
- mc_data_cancel(&thread_ctx->absorbed_by_atmosphere, irealisation);
- mc_data_cancel(&thread_ctx->absorbed_by_receivers, irealisation);
- mc_data_cancel(&thread_ctx->other_absorbed, irealisation);
- mc_data_cancel(&thread_ctx->missing, irealisation);
- mc_data_cancel(&thread_ctx->shadowed, irealisation);
-
- /* Cancel receiver MC estimations */
- htable_receiver_begin(&thread_ctx->mc_rcvs, &r_it);
- htable_receiver_end(&thread_ctx->mc_rcvs, &r_end);
- while (!htable_receiver_iterator_eq(&r_it, &r_end)) {
- struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
- const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
-
- htable_receiver_iterator_next(&r_it);
-
- if (inst->receiver_mask & (int)SSOL_FRONT) {
- cancel_mc_receiver_1side(&mc_rcv->front, irealisation);
- }
- if (inst->receiver_mask & (int)SSOL_BACK) {
- cancel_mc_receiver_1side(&mc_rcv->back, irealisation);
- }
- }
- /* Cancel sampled instance MC estimations */
- htable_sampled_begin(&thread_ctx->mc_samps, &s_it);
- htable_sampled_end(&thread_ctx->mc_samps, &s_end);
- while (!htable_sampled_iterator_eq(&s_it, &s_end)) {
- struct mc_sampled* mc_samp = htable_sampled_iterator_data_get(&s_it);
- htable_sampled_iterator_next(&s_it);
-
- mc_data_cancel(&mc_samp->cos_factor, irealisation);
- mc_data_cancel(&mc_samp->shadowed, irealisation);
-
- /* dst->by_receiver += src->by_receiver; */
- htable_receiver_begin(&mc_samp->mc_rcvs, &r_it);
- htable_receiver_end(&mc_samp->mc_rcvs, &r_end);
- while (!htable_receiver_iterator_eq(&r_it, &r_end)) {
- struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
- const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
- htable_receiver_iterator_next(&r_it);
-
- if (inst->receiver_mask & (int)SSOL_FRONT) {
- cancel_mc_receiver_1side(&mc_rcv->front, irealisation);
- }
- if (inst->receiver_mask & (int)SSOL_BACK) {
- cancel_mc_receiver_1side(&mc_rcv->back, irealisation);
- }
- }
- }
-}
-
/*******************************************************************************
* Exported functions
******************************************************************************/
