solstice-solver

Solver library of the solstice app
git clone git://git.meso-star.com/solstice-solver.git
Log | Files | Refs | README | LICENSE

commit 81152d69d2ecfd82a3e14654f85e416fe10d9d0a
parent fcee059601c6901d5f4592ae3da3e20b9c18df79
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 14 Feb 2017 18:30:09 +0100

First commit, just to save work.

Diffstat:
Msrc/ssol.h | 8++++++++
Msrc/ssol_estimator.c | 142+++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------
Msrc/ssol_estimator_c.h | 154+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----
Msrc/ssol_solver.c | 153+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------
4 files changed, 364 insertions(+), 93 deletions(-)

diff --git a/src/ssol.h b/src/ssol.h @@ -779,6 +779,14 @@ ssol_estimator_get_receiver_status struct ssol_estimator_status* status); SSOL_API res_T +ssol_estimator_get_primary_entity_x_receiver_status + (struct ssol_estimator* estimator, + const struct ssol_instance* prim_instance, + const struct ssol_instance* recv_instance, + const enum ssol_side_flag side, + struct ssol_estimator_status* status); + +SSOL_API res_T ssol_estimator_get_count (const struct ssol_estimator* estimator, size_t* count); diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c @@ -30,7 +30,7 @@ * Helper functions ******************************************************************************/ static res_T -create_per_receiver_mc_data +create_per_prim_recv_mc_data (struct ssol_estimator* estimator, struct ssol_scene* scene) { @@ -45,16 +45,24 @@ create_per_receiver_mc_data const struct ssol_instance* inst = *htable_instance_iterator_data_get(&it); htable_instance_iterator_next(&it); - if(!inst->receiver_mask) continue; + if (inst->receiver_mask) { + res = htable_receiver_set + (&estimator->global_receivers, &inst, &MC_RECV_DATA_NULL); + if (res != RES_OK) goto error; + } - res = htable_receiver_set - (&estimator->global_receivers, &inst, &MC_RECV_DATA_NULL); - if(res != RES_OK) goto error; + if (inst->sample) { + struct mc_per_primary_data data; + init_mc_per_prim_data(estimator->global_primaries.allocator, &data); + res = htable_primary_set(&estimator->global_primaries, &inst, &data); + if (res != RES_OK) goto error; + } } exit: return res; error: htable_receiver_clear(&estimator->global_receivers); + htable_primary_clear(&estimator->global_primaries); goto exit; } @@ -67,6 +75,7 @@ estimator_release(ref_T* ref) ASSERT(ref); dev = estimator->dev; htable_receiver_release(&estimator->global_receivers); + htable_primary_release(&estimator->global_primaries); ASSERT(dev && dev->allocator); MEM_RM(dev->allocator, estimator); SSOL(device_ref_put(dev)); @@ -77,7 +86,7 @@ estimator_release(ref_T* ref) ******************************************************************************/ res_T ssol_estimator_ref_get -(struct ssol_estimator* estimator) + (struct ssol_estimator* estimator) { if (!estimator) return RES_BAD_ARG; ref_get(&estimator->ref); @@ -85,7 +94,7 @@ ssol_estimator_ref_get } res_T -ssol_estimator_ref_put + ssol_estimator_ref_put (struct ssol_estimator* estimator) { if (!estimator) return RES_BAD_ARG; @@ -104,8 +113,8 @@ ssol_estimator_get_status return RES_BAD_ARG; switch (type) { - case SSOL_STATUS_SHADOW: data = &estimator->shadow; break; - case SSOL_STATUS_MISSING: data = &estimator->missing; break; + case SSOL_STATUS_SHADOW: data = &estimator->global_shadow; break; + case SSOL_STATUS_MISSING: data = &estimator->global_missing; break; default: FATAL("Unreachable code.\n"); break; } status->N = estimator->realisation_count; @@ -117,15 +126,13 @@ ssol_estimator_get_status status->irradiance.SE = (status->irradiance.V > 0) ? sqrt(status->irradiance.V / (double)status->N) : 0; - status->absorptivity_loss.E = 0; - status->absorptivity_loss.V = 0; - status->absorptivity_loss.SE = 0; - status->reflectivity_loss.E = 0; - status->reflectivity_loss.V = 0; - status->reflectivity_loss.SE = 0; - status->cos_loss.E = 0; - status->cos_loss.V = 0; - status->cos_loss.SE = 0; + + #define CLEAR_MC_FIELD(Field) \ + status->Field.E = status->Field.V = status->Field.SE = 0 + CLEAR_MC_FIELD(absorptivity_loss); + CLEAR_MC_FIELD(reflectivity_loss); + CLEAR_MC_FIELD(cos_loss); + #undef CLEAR_MC_FIELD return RES_OK; } @@ -148,34 +155,74 @@ ssol_estimator_get_receiver_status status->N = estimator->realisation_count; status->Nf = estimator->failed_count; - status->irradiance.E = data->irradiance.weight / (double)status->N; - status->irradiance.V - = data->irradiance.sqr_weight / (double)status->N - - status->irradiance.E * status->irradiance.E; - status->irradiance.SE - = (status->irradiance.V > 0) ? - sqrt(status->irradiance.V / (double)status->N) : 0; - status->absorptivity_loss.E = data->absorptivity_loss.weight / (double) status->N; - status->absorptivity_loss.V - = data->absorptivity_loss.sqr_weight / (double) status->N - - status->absorptivity_loss.E * status->absorptivity_loss.E; - status->absorptivity_loss.SE - = (status->absorptivity_loss.V > 0) ? - sqrt(status->absorptivity_loss.V / (double) status->N) : 0; - status->reflectivity_loss.E = data->reflectivity_loss.weight / (double) status->N; - status->reflectivity_loss.V - = data->reflectivity_loss.sqr_weight / (double) status->N - - status->reflectivity_loss.E * status->reflectivity_loss.E; - status->reflectivity_loss.SE - = (status->reflectivity_loss.V > 0) ? - sqrt(status->reflectivity_loss.V / (double) status->N) : 0; - status->cos_loss.E = data->cos_loss.weight / (double) status->N; - status->cos_loss.V - = data->cos_loss.sqr_weight / (double) status->N - - status->cos_loss.E * status->cos_loss.E; - status->cos_loss.SE - = (status->cos_loss.V > 0) ? - sqrt(status->cos_loss.V / (double) status->N) : 0; + + #define SET_MC_FIELD(Field) {\ + status->Field.E = data->Field.weight / (double)status->N;\ + status->Field.V\ + = data->Field.sqr_weight / (double) status->N\ + - status->Field.E * status->Field.E;\ + status->Field.SE\ + = (status->Field.V > 0) ?\ + sqrt(status->Field.V / (double) status->N) : 0;\ + } (void) 0 + SET_MC_FIELD(irradiance); + SET_MC_FIELD(absorptivity_loss); + SET_MC_FIELD(reflectivity_loss); + SET_MC_FIELD(cos_loss); + #undef SET_MC_FIELD + + return RES_OK; +} + +res_T +ssol_estimator_get_primary_entity_x_receiver_status + (struct ssol_estimator* estimator, + const struct ssol_instance* prim_instance, + const struct ssol_instance* recv_instance, + const enum ssol_side_flag side, + struct ssol_estimator_status* status) +{ + struct mc_per_primary_data* prim_data = NULL; + struct mc_per_receiver_1side_data* data = NULL; + if (!estimator || !prim_instance || !recv_instance || !status + || (side != SSOL_BACK && side != SSOL_FRONT)) + return RES_BAD_ARG; + + /* Check if prim_instance is a primary entity */ + prim_data = estimator_get_primary_entity_data + (&estimator->global_primaries, prim_instance); + if (prim_data == NULL) return RES_BAD_ARG; + + /* Check if a receiver is defined for this instance/side */ + data = estimator_get_prim_recv_data(prim_data, recv_instance, side); + if (data == NULL) return RES_BAD_ARG; + + /* realisation count for this primary */ + status->N = prim_data->nb_samples; + status->Nf = prim_data->nb_failed; + + #define CLEAR_MC_FIELD(Field) \ + status->Field.E = status->Field.V = status->Field.SE = 0 + #define SET_MC_FIELD(Field) {\ + if (status->N > 0) {\ + status->Field.E = data->Field.weight / (double)status->N;\ + status->Field.V\ + = data->Field.sqr_weight / (double) status->N\ + - status->Field.E * status->Field.E;\ + status->Field.SE\ + = (status->Field.V > 0) ?\ + sqrt(status->Field.V / (double) status->N) : 0;\ + } else {\ + CLEAR_MC_FIELD(Field);\ + }\ + } (void) 0 + SET_MC_FIELD(irradiance); + SET_MC_FIELD(absorptivity_loss); + SET_MC_FIELD(reflectivity_loss); + SET_MC_FIELD(cos_loss); + #undef CLEAR_MC_FIELD + #undef SET_MC_FIELD + return RES_OK; } @@ -241,11 +288,12 @@ estimator_create } htable_receiver_init(dev->allocator, &estimator->global_receivers); + htable_primary_init(dev->allocator, &estimator->global_primaries); SSOL(device_ref_get(dev)); estimator->dev = dev; ref_init(&estimator->ref); - res = create_per_receiver_mc_data(estimator, scene); + res = create_per_prim_recv_mc_data(estimator, scene); if(res != RES_OK) goto error; exit: diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h @@ -20,6 +20,11 @@ #include <rsys/hash_table.h> struct mem_allocator; +struct ssol_instance; + +/* + * Common stuff + */ static FINLINE int side_idx(const enum ssol_side_flag side) @@ -37,6 +42,10 @@ struct mc_data { #define MC_DATA_NULL__ { 0, 0 } static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__; +/* + * Stuff for by-receiver results + */ + struct mc_per_receiver_1side_data { struct mc_data irradiance; struct mc_data absorptivity_loss; @@ -55,7 +64,8 @@ struct mc_per_receiver_data { struct mc_per_receiver_1side_data back; }; -#define MC_RECV_DATA_NULL__ { MC_RECV_1SIDE_DATA_NULL__, MC_RECV_1SIDE_DATA_NULL__ } +#define MC_RECV_DATA_NULL__ {\ + MC_RECV_1SIDE_DATA_NULL__, MC_RECV_1SIDE_DATA_NULL__ } static const struct mc_per_receiver_data MC_RECV_DATA_NULL = MC_RECV_DATA_NULL__; @@ -71,22 +81,130 @@ init_mc_per_recv_data } /* Define the htable_receiver data structure */ -struct ssol_instance; #define HTABLE_NAME receiver #define HTABLE_KEY const struct ssol_instance* #define HTABLE_DATA struct mc_per_receiver_data #define HTABLE_FUNCTOR_INIT init_mc_per_recv_data #include <rsys/hash_table.h> +#undef HTABLE_NAME +#undef HTABLE_KEY +#undef HTABLE_DATA +#undef HTABLE_FUNCTOR_INIT + +/* + * Stuff for by-primary-entity results + */ + +struct mc_per_primary_data { + /* global data for this entity */ + struct mc_data cos_loss; + struct mc_data shadow_loss; + double area; + double base_sun_cos; + size_t nb_samples; + size_t nb_failed; + /* by-receptor data for this entity */ + struct htable_receiver by_receiver; +}; + +static INLINE void +init_mc_per_prim_data + (struct mem_allocator* alloc, + struct mc_per_primary_data* data) +{ + ASSERT(alloc && data); + data->area = 0; + data->base_sun_cos = 0; + data->cos_loss = MC_DATA_NULL; + data->shadow_loss = MC_DATA_NULL; + data->nb_samples = 0; + data->nb_failed = 0; + htable_receiver_init(alloc, &data->by_receiver); +} + +static INLINE void +release_mc_per_prim_data(struct mc_per_primary_data* data) +{ + ASSERT(data); + htable_receiver_release(&data->by_receiver); +} + +#define PRIM_COPY(Dst, Src) {\ + Dst->area = Src->area;\ + Dst->base_sun_cos = Src->base_sun_cos;\ + Dst->nb_failed = Src->nb_failed;\ + Dst->nb_samples = Src->nb_samples;\ + Dst->shadow_loss = Src->shadow_loss;\ + Dst->cos_loss = Src->cos_loss;\ +} (void)0 + +static INLINE res_T +copy_mc_per_prim_data + (struct mc_per_primary_data* dst, + const struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy(&dst->by_receiver, &src->by_receiver); +} + +static INLINE res_T +copy_and_release_mc_per_prim_data + (struct mc_per_primary_data* dst, + struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy_and_release(&dst->by_receiver, &src->by_receiver); +} + +static INLINE res_T +copy_and_clear_mc_per_prim_data + (struct mc_per_primary_data* dst, + struct mc_per_primary_data* src) +{ + ASSERT(dst && src); + PRIM_COPY(dst, src); + return htable_receiver_copy_and_clear(&dst->by_receiver, &src->by_receiver); +} + +#undef PRIM_COPY + +/* Define the htable_primary data structure */ +#define HTABLE_NAME primary +#define HTABLE_KEY const struct ssol_instance* +#define HTABLE_DATA struct mc_per_primary_data +#define HTABLE_DATA_FUNCTOR_INIT init_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_RELEASE release_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY copy_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE copy_and_release_mc_per_prim_data +#define HTABLE_DATA_FUNCTOR_COPY_AND_CLEAR copy_and_clear_mc_per_prim_data +#include <rsys/hash_table.h> +#undef HTABLE_NAME +#undef HTABLE_KEY +#undef HTABLE_DATA +#undef HTABLE_DATA_FUNCTOR_INIT +#undef HTABLE_DATA_FUNCTOR_RELEASE +#undef HTABLE_DATA_FUNCTOR_COPY +#undef HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE +#undef HTABLE_DATA_FUNCTOR_COPY_AND_CLEAR + +/* + * ssol_estimator + * gathers results by receiver, by primary entity, and by primary-receiver couple + */ struct ssol_estimator { size_t realisation_count; size_t failed_count; - /* the implicit MC computations */ - struct mc_data shadow; - struct mc_data missing; - /* 1 global MC per receiver */ + /* implicit, global MC computations */ + struct mc_data global_shadow; + struct mc_data global_missing; + /* per-receiver MC computations */ struct htable_receiver global_receivers; - /* areas */ + /* per-primary entity MC computations */ + struct htable_primary global_primaries; + /* scene areas */ double sampled_area, primary_area; struct ssol_device* dev; @@ -113,4 +231,26 @@ estimator_get_receiver_data return side == SSOL_FRONT ? &data->front : &data->back; } +static FINLINE struct mc_per_primary_data* +estimator_get_primary_entity_data +(struct htable_primary* primaries, + const struct ssol_instance* instance) +{ + struct mc_per_primary_data* data; + ASSERT(primaries && instance); + if (!instance->sample) return NULL; + data = htable_primary_find(primaries, &instance); + return data; +} + +static FINLINE struct mc_per_receiver_1side_data* +estimator_get_prim_recv_data + (struct mc_per_primary_data* primary_data, + const struct ssol_instance* instance, + const enum ssol_side_flag side) +{ + ASSERT(primary_data && instance); + return estimator_get_receiver_data(&primary_data->by_receiver, instance, side); +} + #endif /* SSOL_ESTIMATOR_C_H */ diff --git a/src/ssol_solver.c b/src/ssol_solver.c @@ -47,6 +47,7 @@ struct point { const struct ssol_instance* inst; const struct shaded_shape* sshape; + struct mc_per_primary_data* primary; struct s3d_primitive prim; double N[3]; double pos[3]; @@ -63,6 +64,7 @@ struct point { #define POINT_NULL__ { \ NULL, /* Instance */ \ NULL, /* Shaded shape */ \ + NULL, /* Primary data */ \ S3D_PRIMITIVE_NULL__, /* Primitive */ \ {0, 0, 0}, /* Normal */ \ {0, 0, 0}, /* Position */ \ @@ -82,7 +84,7 @@ static int point_init (struct point* pt, struct ssol_scene* scn, - const double sampled_area, + struct ssol_estimator* estimator, struct s3d_scene_view* view_samp, struct s3d_scene_view* view_rt, struct ranst_sun_dir* ran_sun_dir, @@ -118,8 +120,8 @@ point_init /* Initialise the Monte Carlo weight */ cos_sun = fabs(d3_dot(pt->N, pt->dir)); - pt->weight = scn->sun->dni * sampled_area * cos_sun; - pt->cos_loss = scn->sun->dni * sampled_area * (1 - cos_sun); + pt->weight = scn->sun->dni * estimator->sampled_area * cos_sun; + pt->cos_loss = scn->sun->dni * estimator->sampled_area * (1 - cos_sun); pt->absorptivity_loss = pt->reflectivity_loss = 0; /* Retrieve the sampled instance and shaded shape */ @@ -129,6 +131,21 @@ point_init pt->sshape = darray_shaded_shape_cdata_get (&pt->inst->object->shaded_shapes) + id; + /* store primary entity related weights */ + pt->primary = + estimator_get_primary_entity_data(&estimator->global_primaries, pt->inst); + if (!pt->primary) { + struct mc_per_primary_data data; + init_mc_per_prim_data(estimator->global_primaries.allocator, &data); + htable_primary_set(&estimator->global_primaries, &pt->inst, &data); + pt->primary = + estimator_get_primary_entity_data(&estimator->global_primaries, pt->inst); + ASSERT(pt->primary); + } + pt->primary->cos_loss.weight += pt->cos_loss; + pt->primary->cos_loss.sqr_weight += pt->cos_loss * pt->cos_loss; + pt->primary->nb_samples++; + /* For punched surface, retrieve the sampled position and normal onto the * quadric surface */ if(pt->sshape->shape->type == SHAPE_PUNCHED) { @@ -339,7 +356,8 @@ ssol_solve FILE* output, struct ssol_estimator** out_estimator) { - struct htable_receiver_iterator it, end; + struct htable_receiver_iterator r_it, r_end; + struct htable_primary_iterator p_it, p_end; struct s3d_scene_view* view_rt = NULL; struct s3d_scene_view* view_samp = NULL; struct s3d_scene_view* view_prim = NULL; @@ -350,11 +368,13 @@ ssol_solve struct ssp_rng_proxy* rng_proxy = NULL; struct mc_data* mc_shadows = NULL; struct mc_data* mc_missings = NULL; + size_t* failed_counts = NULL; struct htable_receiver* mc_rcvs = NULL; + struct htable_primary* mc_prims = NULL; struct ssol_estimator* estimator = NULL; float area; int nthreads = 0; - int nrealisations = 0; + int64_t nrealisations = 0; int i = 0; ATOMIC res = RES_OK; ASSERT(nrealisations <= INT_MAX); @@ -367,7 +387,7 @@ ssol_solve /* CL compiler supports OpenMP parallel loop whose indices are signed. The * following line ensures that the unsigned number of realisations does not * overflow the realisation index. */ - if(realisations_count > INT_MAX) { + if(realisations_count > INT64_MAX) { res = RES_BAD_ARG; goto error; } @@ -410,6 +430,8 @@ ssol_solve CREATE(mc_shadows); CREATE(mc_missings); CREATE(mc_rcvs); + CREATE(mc_prims); + CREATE(failed_counts); #undef CREATE /* Setup per thread data structures */ @@ -421,12 +443,24 @@ ssol_solve htable_receiver_init(scn->dev->allocator, mc_rcvs + i); res = htable_receiver_copy(mc_rcvs + i, &estimator->global_receivers); if(res != RES_OK) goto error; - htable_receiver_begin(mc_rcvs + i, &it); - htable_receiver_end(mc_rcvs + i, &end); - while (!htable_receiver_iterator_eq(&it, &end)) { - struct mc_per_receiver_data* estimator_data = htable_receiver_iterator_data_get(&it); + htable_receiver_begin(mc_rcvs + i, &r_it); + htable_receiver_end(mc_rcvs + i, &r_end); + while (!htable_receiver_iterator_eq(&r_it, &r_end)) { + struct mc_per_receiver_data* estimator_data + = htable_receiver_iterator_data_get(&r_it); *estimator_data = MC_RECV_DATA_NULL; - htable_receiver_iterator_next(&it); + htable_receiver_iterator_next(&r_it); + } + htable_primary_init(scn->dev->allocator, mc_prims + i); + res = htable_primary_copy(mc_prims + i, &estimator->global_primaries); + if (res != RES_OK) goto error; + htable_primary_begin(mc_prims + i, &p_it); + htable_primary_end(mc_prims + i, &p_end); + while (!htable_primary_iterator_eq(&p_it, &p_end)) { + struct mc_per_primary_data* estimator_data + = htable_primary_iterator_data_get(&p_it); + init_mc_per_prim_data(scn->dev->allocator, estimator_data); + htable_primary_iterator_next(&p_it); } } @@ -455,10 +489,12 @@ ssol_solve receiver = mc_rcvs + ithread; /* Find a new starting point of the radiative random walk */ - is_lit = point_init(&pt, scn, estimator->sampled_area, view_samp, view_rt, + is_lit = point_init(&pt, scn, estimator, view_samp, view_rt, ran_sun_dir, ran_sun_wl, rng); if(!is_lit) { /* The starting point is not lit */ + pt.primary->shadow_loss.weight += pt.weight; + pt.primary->shadow_loss.sqr_weight += pt.weight; shadow->weight += pt.weight; shadow->sqr_weight += pt.weight * pt.weight; continue; @@ -483,14 +519,17 @@ ssol_solve } data = estimator_get_receiver_data(receiver, pt.inst, pt.side); ASSERT(data); - data->irradiance.weight += pt.weight; - data->irradiance.sqr_weight += pt.weight*pt.weight; - data->absorptivity_loss.weight += pt.absorptivity_loss; - data->absorptivity_loss.sqr_weight += pt.absorptivity_loss * pt.absorptivity_loss; - data->reflectivity_loss.weight += pt.reflectivity_loss; - data->reflectivity_loss.sqr_weight += pt.reflectivity_loss * pt.reflectivity_loss; - data->cos_loss.weight += pt.cos_loss; - data->cos_loss.sqr_weight += pt.cos_loss * pt.cos_loss; + #define MC_ADD(Field) MC_ADD2(Field, Field) + #define MC_ADD2(Data_Field, Pt_Field) {\ + data->Data_Field.weight += pt.Pt_Field;\ + data->Data_Field.sqr_weight += pt.Pt_Field * pt.Pt_Field;\ + } (void)0 + MC_ADD2(irradiance, weight); + MC_ADD(absorptivity_loss); + MC_ADD(reflectivity_loss); + MC_ADD(cos_loss); + #undef MC_ADD + #undef MC_ADD2 hit_a_receiver = 1; } @@ -549,31 +588,36 @@ ssol_solve } } - /* Merge per thread estimations */ + /* + * Merge per thread estimations + */ + + /* global results */ FOR_EACH(i, 0, nthreads) { - estimator->shadow.weight += mc_shadows[i].weight; - estimator->shadow.sqr_weight += mc_shadows[i].sqr_weight; - estimator->missing.weight += mc_missings[i].weight; - estimator->missing.sqr_weight += mc_missings[i].sqr_weight; + estimator->global_shadow.weight += mc_shadows[i].weight; + estimator->global_shadow.sqr_weight += mc_shadows[i].sqr_weight; + estimator->global_missing.weight += mc_missings[i].weight; + estimator->global_missing.sqr_weight += mc_missings[i].sqr_weight; + estimator->failed_count += failed_counts[i]; } - /* Merge per thread estimations */ - htable_receiver_begin(&estimator->global_receivers, &it); - htable_receiver_end(&estimator->global_receivers, &end); - while (!htable_receiver_iterator_eq(&it, &end)) { - struct mc_per_receiver_data* estimator_data - = htable_receiver_iterator_data_get(&it); - const struct ssol_instance* key = *htable_receiver_iterator_key_get(&it); - htable_receiver_iterator_next(&it); + /* per-receiver results */ + htable_receiver_begin(&estimator->global_receivers, &r_it); + htable_receiver_end(&estimator->global_receivers, &r_end); + while (!htable_receiver_iterator_eq(&r_it, &r_end)) { + struct mc_per_receiver_data* est_data + = htable_receiver_iterator_data_get(&r_it); + const struct ssol_instance* key = *htable_receiver_iterator_key_get(&r_it); + htable_receiver_iterator_next(&r_it); FOR_EACH(i, 0, nthreads) { /* Sum both sides, even if no receiver is defined to avoid tests */ - struct mc_per_receiver_data* thread_data + struct mc_per_receiver_data* th_data = htable_receiver_find(mc_rcvs + i, &key); - #define ACCUM_WEIGHT(Name) { \ - estimator_data->front.Name.weight += thread_data->front.Name.weight; \ - estimator_data->back.Name.weight += thread_data->front.Name.weight; \ - estimator_data->front.Name.sqr_weight += thread_data->front.Name.sqr_weight;\ - estimator_data->back.Name.sqr_weight += thread_data->front.Name.sqr_weight;\ + #define ACCUM_WEIGHT(Field) { \ + est_data->front.Field.weight += th_data->front.Field.weight; \ + est_data->back.Field.weight += th_data->back.Field.weight; \ + est_data->front.Field.sqr_weight += th_data->front.Field.sqr_weight;\ + est_data->back.Field.sqr_weight += th_data->back.Field.sqr_weight; \ } (void)0 ACCUM_WEIGHT(irradiance); ACCUM_WEIGHT(absorptivity_loss); @@ -582,6 +626,32 @@ ssol_solve #undef ACCUM_WEIGHT } } + + /* per-primary entity results */ + htable_primary_begin(&estimator->global_primaries, &p_it); + htable_primary_end(&estimator->global_primaries, &p_end); + while (!htable_primary_iterator_eq(&p_it, &p_end)) { + struct mc_per_primary_data* pri_data + = htable_primary_iterator_data_get(&p_it); + const struct ssol_instance* key = *htable_primary_iterator_key_get(&p_it); + htable_primary_iterator_next(&p_it); + FOR_EACH(i, 0, nthreads) { + /* Sum both sides, even if no receiver is defined to avoid tests */ + struct mc_per_primary_data* th_data + = htable_primary_find(mc_prims + i, &key); + #define ACCUM_MC_WEIGHT(Field) { \ + pri_data->Field.weight += th_data->Field.weight; \ + pri_data->Field.sqr_weight += th_data->Field.sqr_weight; \ + } (void)0 + ACCUM_MC_WEIGHT(cos_loss); + ACCUM_MC_WEIGHT(shadow_loss); + pri_data->nb_samples += th_data->nb_samples; + pri_data->nb_failed += th_data->nb_failed; + //pri_data->by_receiver += th_data->by_receiver; + #undef ACCUM_MC_WEIGHT + } + } + estimator->realisation_count += realisations_count; exit: @@ -603,8 +673,13 @@ exit: FOR_EACH(i, 0, nthreads) htable_receiver_release(mc_rcvs + i); sa_release(mc_rcvs); } + if (mc_prims) { + FOR_EACH(i, 0, nthreads) htable_primary_release(mc_prims + i); + sa_release(mc_prims); + } sa_release(mc_shadows); sa_release(mc_missings); + sa_release(failed_counts); if(out_estimator) *out_estimator = estimator; return (res_T)res; error: