commit 226cb3f14e393e9186ecb67d85b20a6e6e0d9035
parent 27df4af5d16d6ed2fd96afbaa246d94e58a95c1a
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Wed, 15 Oct 2025 16:17:06 +0200
Setup darray_accum to handle weight accumulation in MVREA computation
This commit refactors how sphor_compute_MVREA handles weight
accumulators. The introduction of the struct accum and associated data
structures enables the management of multiple random variables within
the same algorithm. The list of all possible accumulators can now be
easily created and initialized by iterating over all sensors in the
scene.
A separate copy of this list is created for each thread involved in the
computation, removing any responsibility for multithreading from the
realization function. The aggregation of per-thread accumulators into
the global accumulator is performed at the end of the realization loop,
in the same spirit of what is usually done when the weight is a simple
scalar (double).
This commit also extends the internal API of the accum module by
extracting logic that may be reused in different contexts, such as:
1. Aggregating per-thread accumulators into the main accumulator
(accum_threads_aggregate)
2. Registering new accumulators with custom string fields
(register_accum)
3. Updating accumulator values dynamically using an on-the-fly
constructed accum_key (update_accum)
Diffstat:
5 files changed, 312 insertions(+), 48 deletions(-)
diff --git a/src/sphor.c b/src/sphor.c
@@ -148,13 +148,11 @@ sphor_ref_put(struct sphor* sphor)
res_T
sphor_run(struct sphor* sphor)
{
- double estim = 0;
- double std = 0;
res_T res = RES_OK;
if (NULL == sphor) { goto error; }
- res = sphor_compute_MVREA(sphor, &estim, &std);
+ res = sphor_compute_MVREA(sphor);
if (RES_OK != res) { goto error; }
exit:
diff --git a/src/sphor_accum.c b/src/sphor_accum.c
@@ -94,3 +94,75 @@ accum_copy_and_release
res = str_copy_and_release(&dst->component, &src->component);
return res;
}
+
+extern LOCAL_SYM res_T
+register_accum
+ (struct mem_allocator* allocator,
+ char* observable,
+ char* sensor,
+ char* component,
+ struct darray_accum *accums)
+{
+ struct accum accum = ACCUM_NULL;
+ res_T res = RES_OK;
+
+ accum_init(allocator, &accum);
+
+ str_set(&accum.observable, observable);
+ str_set(&accum.sensor, sensor);
+ str_set(&accum.component, component);
+
+ res = darray_accum_push_back(accums, &accum);
+ if (res != RES_OK) { goto error; }
+exit:
+ accum_release(&accum);
+ return res;
+error:
+ goto exit;
+}
+
+extern LOCAL_SYM void
+accum_threads_aggregate
+ (size_t nthreads,
+ struct darray_accum* accums_threads,
+ struct darray_accum* accums)
+{
+ size_t i = 0, j = 0;
+
+ ASSERT(NULL != accums_threads);
+ ASSERT(NULL != accums);
+
+ FOR_EACH(i, 0, darray_accum_size_get(accums)){
+ struct accum *main_accum;
+ main_accum = &darray_accum_data_get(accums)[i];
+ FOR_EACH(j, 0, nthreads) {
+ struct accum thread_accum;
+ thread_accum = darray_accum_data_get(&accums_threads[j])[i];
+ main_accum->sum += thread_accum.sum;
+ main_accum->sum2 += thread_accum.sum2;
+ }
+ }
+}
+
+extern LOCAL_SYM void
+update_accum
+ (const struct htable_accum2id* accum2id,
+ const struct accum_key* accum_key,
+ const double weight,
+ struct darray_accum* accums)
+{
+ size_t* paccum_id = NULL;
+ struct accum* accum = NULL;
+
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+
+
+ paccum_id = htable_accum2id_find((struct htable_accum2id *)accum2id, accum_key);
+ ASSERT(NULL != paccum_id);
+
+ accum = &darray_accum_data_get(accums)[*paccum_id];
+
+ accum->sum += weight;
+ accum->sum2 += weight*weight;
+}
diff --git a/src/sphor_accum.h b/src/sphor_accum.h
@@ -32,8 +32,8 @@
struct accum {
struct str observable; /* Type of observable */
struct str sensor; /* Volume name or surface name */
- struct str component; /* Optional: used to discriminate part of the events (ex:
- prop_rad name in volumes) */
+ struct str component; /* Optional: used to discriminate part of the events
+ (ex: prop_rad name in volumes) */
double sum;
double sum2;
size_t n_realizations;
@@ -134,4 +134,25 @@ accum_key_hash
/* Include rsys/hash_table a second time to generate the hash table */
#include <rsys/hash_table.h>
+extern LOCAL_SYM res_T
+register_accum
+ (struct mem_allocator* allocator,
+ char* observable,
+ char* sensor,
+ char* component,
+ struct darray_accum *accums);
+
+extern LOCAL_SYM void
+accum_threads_aggregate
+ (size_t nthreads,
+ struct darray_accum* accums_threads,
+ struct darray_accum* accums);
+
+extern LOCAL_SYM void
+update_accum
+ (const struct htable_accum2id* accum2id,
+ const struct accum_key* accum_key,
+ const double weight,
+ struct darray_accum* accums);
+
#endif /* SPHOR_ACCUM_H */
diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c
@@ -24,6 +24,7 @@
#include "sphor.h"
#include "sphor_c.h"
+#include "sphor_accum.h"
#include "sphor_compute_mvrea.h"
#include "sphor_interface.h"
#include "sphor_ran_source.h"
@@ -32,6 +33,8 @@
#include <star/sphin.h>
#include <star/s3d.h>
#include <star/ssp.h>
+#include <rsys/dynamic_array.h>
+#include <rsys/hash_table.h>
#include <omp.h>
enum ray_interface_interaction_type {
@@ -45,6 +48,147 @@ enum ray_interface_interaction_type {
* Helper functions
******************************************************************************/
static res_T
+setup_MVREA_accums_surfaces
+ (struct sphor* sphor,
+ struct htable_accum2id *accum2id,
+ struct darray_accum* accums)
+{
+ char* surface_name = NULL;
+ size_t accum_id = 0;
+ size_t surface_count = 0;
+ size_t i_surface = 0;
+ struct sphin_config* config = sphor->config;
+ struct sphin_surface* surface = NULL;
+ struct sphin_sensor_surface* sensor = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+
+ res = sphin_config_get_surface_count(config, &surface_count);
+ if (RES_OK != res){ goto error; }
+
+ FOR_EACH(i_surface, 0, surface_count){
+ res = sphin_config_get_surface(config, i_surface, &surface);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_surface_get_sensor(surface, &sensor);
+ if (RES_OK != res){ goto error; }
+ if (NULL != sensor) {
+ res = sphin_surface_get_name(surface, &surface_name);
+ if (RES_OK != res){ goto error; }
+
+ res = register_accum
+ (sphor->allocator, "LOSSES", surface_name, NULL, accums);
+ if (RES_OK != res){ goto error; }
+ accum_id = darray_accum_size_get(accums) - 1;
+
+ /* Set corresponding entry in the accum2id hash table */
+ struct accum_key accum_key = ACCUM_KEY_NULL;
+ accum_key.sensor_type = SPHOR_SENSOR_SURFACE;
+ accum_key.sensor_id = i_surface;
+ res = htable_accum2id_set(accum2id, &accum_key, &accum_id);
+ if (RES_OK != res){ goto error; }
+ }
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_MVREA_accums_volumes
+ (struct sphor* sphor,
+ struct htable_accum2id *accum2id,
+ struct darray_accum* accums)
+{
+ char* volume_name = NULL;
+ char* prop_rad_name = NULL;
+ size_t accum_id = 0;
+ size_t volume_count = 0;
+ size_t prop_rad_count = 0;
+ size_t i_volume = 0, j_proprad = 0;
+ struct sphin_config* config = sphor->config;
+ struct sphin_volume* volume = NULL;
+ struct sphin_sensor_volume* sensor = NULL;
+ struct sphin_prop_rad* prop_rad = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+
+ res = sphin_config_get_volume_count(config, &volume_count);
+ if (RES_OK != res){ goto error; }
+
+ FOR_EACH(i_volume, 0, volume_count){
+ res = sphin_config_get_volume(config, i_volume, &volume);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_volume_get_sensor(volume, &sensor);
+ if (RES_OK != res){ goto error; }
+ if (NULL != sensor) {
+ res = sphin_volume_get_name(volume, &volume_name);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_volume_get_prop_rad_count(volume, &prop_rad_count);
+ if (RES_OK != res){ goto error; }
+
+ FOR_EACH(j_proprad, 0, prop_rad_count){
+
+ res = sphin_volume_get_prop_rad(volume, j_proprad, &prop_rad);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_prop_rad_get_name(prop_rad, &prop_rad_name);
+ if (RES_OK != res){ goto error; }
+
+ res = register_accum
+ (sphor->allocator, "MVREA", volume_name, prop_rad_name, accums);
+ if (RES_OK != res){ goto error; }
+ accum_id = darray_accum_size_get(accums) - 1;
+
+ /* Set corresponding entry in the accum2id hash table */
+ struct accum_key accum_key = ACCUM_KEY_NULL;
+ accum_key.sensor_type = SPHOR_SENSOR_VOLUME;
+ accum_key.sensor_id = i_volume;
+ accum_key.component_id = j_proprad;
+ res = htable_accum2id_set(accum2id, &accum_key, &accum_id);
+ if (RES_OK != res){ goto error; }
+ }
+ }
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_MVREA_accums
+ (struct sphor *sphor,
+ struct htable_accum2id *accum2id,
+ struct darray_accum *accums)
+{
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accums);
+
+ res = setup_MVREA_accums_volumes(sphor, accum2id, accums);
+ if (RES_OK != res){ goto error; }
+
+ res = setup_MVREA_accums_surfaces(sphor, accum2id, accums);
+ if (RES_OK != res){ goto error; }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
volume_compute_total_ka
(struct sphor* sphor,
struct interface interface,
@@ -256,8 +400,8 @@ static res_T
compute_MVREA_realization
(struct sphor* sphor,
struct ssp_rng* rng,
- double* weight/* TODO make this a list in sphor, since we can have multiple
- sensors. Each sensor should have its own weight */)
+ struct htable_accum2id* accum2id,
+ struct darray_accum* accums)
{
double dir[3] = {0}; /* Current direction in path sampling */
double pos[3] = {0}; /* Current position in path sampling */
@@ -268,6 +412,7 @@ compute_MVREA_realization
enum ray_interface_interaction_type ray_interface_interaction_type =
RAY_INTERFACE_INTERACTION_NONE__;
enum sphin_side side_id = SPHIN_SIDE_NONE__;
+ struct accum_key accum_key = ACCUM_KEY_NULL;
struct interface interface = INTERFACE_NULL;
struct s3d_hit hit = S3D_HIT_NULL;
struct s3d_primitive prim = S3D_PRIMITIVE_NULL;
@@ -302,7 +447,7 @@ compute_MVREA_realization
/* If there is no intersection in the sampled direction from the sampled
* position, the photon is lost and it counts zero in the weight. */
- if (S3D_HIT_NONE(&hit)) { *weight = 0; break; }
+ if (S3D_HIT_NONE(&hit)) { break; }
/* Sample the free path length to absorption from an exponential
* distribution with rate parameter ka */
@@ -329,11 +474,15 @@ compute_MVREA_realization
res = sphin_sensor_volume_get_response_function
(sensor_volume, &response_function);
if (RES_OK != res){ goto error; }
- *weight = response_function;
+
+ /* Update accumulator corresponding to the photon absortion */
+ accum_key.sensor_type = SPHOR_SENSOR_VOLUME;
+ accum_key.sensor_id = interface.volumes[side_id];
+ accum_key.component_id = 0; /* TODO HARD CODED TO TEST: CHANGE IT */
+ update_accum(accum2id, &accum_key, response_function, accums);
}
/* If the volume is not a sensor, the absorbed photon does not count in
- * the weight */
- else{ *weight = 0; }
+ * the weight. The weight is implictly 0 */
break;
}
@@ -349,7 +498,7 @@ compute_MVREA_realization
switch (ray_interface_interaction_type) {
/* An absorption took place, the weight is zero and the path stops */
case RAY_INTERFACE_INTERACTION_ABSORPTION:
- *weight = 0;
+ /* TODO: verify if surface is sensor and update its accum */
stop = 1;
break;
/* Else, if a reflection takes place, sample a new direction according to
@@ -390,17 +539,16 @@ error:
******************************************************************************/
res_T
sphor_compute_MVREA
- (struct sphor* sphor,
- double* out_estim,
- double* out_std)
+ (struct sphor* sphor)
{
res_T res = RES_OK;
size_t nthreads, samples;
size_t i; /* iterator */
struct ssp_rng_proxy *rng_proxy = NULL;
struct ssp_rng **rngs = NULL;
- double* sum_thread = NULL;
- double* sum2_thread = NULL;
+ struct darray_accum accums;
+ struct htable_accum2id accum2id;
+ struct darray_accum *accums_threads = NULL;
double sum = 0, sum2 = 0;
double estim = 0, std = 0;
@@ -410,67 +558,94 @@ sphor_compute_MVREA
nthreads = (size_t)omp_get_num_procs();
}
- /* Création du générateur mandataire RNG_MT19937_64 (Mersenne Twister) */
+ /* Initialize the global accumulator and the hash table that allows to
+ * retrieve the id of a given accumulator with an accum_key */
+ darray_accum_init(sphor->allocator, &accums);
+ htable_accum2id_init(sphor->allocator, &accum2id);
+
+ /* Create and initialize dynamic array o accums (one per observable) */
+ res = setup_MVREA_accums(sphor, &accum2id, &accums);
+ if (RES_OK != res){ goto error; }
+
+ /* Create an array containing nthreads pointers to different accums lists */
+ accums_threads = mem_calloc(nthreads, sizeof(*accums_threads));
+ if (NULL == accums_threads) { res = RES_MEM_ERR; goto error; }
+
+ /* Mem allocation of the accummulator list of each one of the threads
+ * Each thread handles a copy of the global accum list just created */
+ FOR_EACH(i, 0, nthreads) {
+ darray_accum_init(sphor->allocator, &accums_threads[i]);
+ res = darray_accum_copy(&accums_threads[i], &accums);
+ if (RES_OK != res){ goto error; }
+ }
+
+ /* Create of the proxy generator RNG_MT19937_64 (Mersenne Twister) */
res = ssp_rng_proxy_create(NULL, SSP_RNG_MT19937_64, nthreads, &rng_proxy);
if (res != RES_OK) { goto error; }
- /* Allocation des generateurs aleatoire pour chaque processus */
+ /* Create an array containing nthreads pointers to different rngs */
rngs = mem_calloc(nthreads, sizeof(*rngs));
- if (NULL == rngs) {
- res = RES_MEM_ERR;
- goto error;
- }
- /* Le generateur mandataire attribue un generateur par sequence aleatoire
- unique que nous stockons dans le tableau rngs (un par processus) */
+ if (NULL == rngs) { res = RES_MEM_ERR; goto error; }
+
+ /* Set one generator per thread using the proxy generator */
FOR_EACH(i, 0, nthreads) {
res = ssp_rng_proxy_create_rng(rng_proxy, i, &rngs[i]);
- if (RES_OK != res) { goto error; }
+ if (RES_OK != res){ goto error; }
}
omp_set_num_threads((int)nthreads);
- sum_thread = mem_calloc(nthreads, sizeof(*sum_thread));
- sum2_thread = mem_calloc(nthreads, sizeof(*sum2_thread));
+
/* Realizations loop */
#pragma omp parallel for schedule(static)
for(i=0; i<samples; i++) {
const int ithread = omp_get_thread_num();
- double w = 0;
res_T res_local = RES_OK;
if (RES_OK != res) continue;
- res_local = compute_MVREA_realization(sphor, rngs[ithread], &w);
+ res_local = compute_MVREA_realization
+ (sphor, rngs[ithread], &accum2id, &accums_threads[ithread]);
if (RES_OK != res_local) {
res = res_local;
- } else {
- sum_thread[ithread] += w;
- sum2_thread[ithread] += w*w;
}
}
- if(res != RES_OK) goto error;
+ if (res != RES_OK) { goto error; }
- FOR_EACH(i, 0, nthreads) {
- sum += sum_thread[i];
- sum2 += sum2_thread[i];
+ /* Aggregate accums from the different threads in the main accumulator */
+ accum_threads_aggregate(nthreads, accums_threads, &accums);
+
+ /* TEMP: aggregate all the different accums sums in sum and sum2 to test the
+ * algorithm. TODO replace by output function */
+ FOR_EACH(i, 0, darray_accum_size_get(&accums)){
+ struct accum* accum;
+ accum = darray_accum_data_get(&accums) + i;
+ sum += accum->sum;
+ sum2 += accum->sum2;
}
estim = sum/(double)samples;
std = sqrt
((sum2/(double)samples - estim*estim)/((double)samples-1));
-
printf("%f +/- %f\n", estim, std);
exit:
- if(sum_thread) mem_rm(sum_thread);
- if(sum2_thread) mem_rm(sum2_thread);
- /* Liberation de l'ensemble des generateurs */
- if(rng_proxy) ssp_rng_proxy_ref_put(rng_proxy);
- if(rngs) {
+ /* Free each one of the thread accumulators */
+ if (accums_threads) {
+ FOR_EACH(i, 0, nthreads) {
+ darray_accum_release(&accums_threads[i]);
+ }
+ mem_rm(accums_threads);
+ }
+ /* Free the proxy generator and each one of the thread generators */
+ if (rng_proxy) ssp_rng_proxy_ref_put(rng_proxy);
+ if (rngs) {
FOR_EACH(i, 0, nthreads) {
if(rngs[i]) ssp_rng_ref_put(rngs[i]);
}
mem_rm(rngs);
}
- *out_estim = estim;
- *out_std = std;
+ /* Free local data structures accums and accum2id */
+ darray_accum_release(&accums);
+ htable_accum2id_release(&accum2id);
+
return res;
error:
goto exit;
diff --git a/src/sphor_compute_mvrea.h b/src/sphor_compute_mvrea.h
@@ -32,8 +32,6 @@ struct sphor_create_args;
extern LOCAL_SYM res_T
sphor_compute_MVREA
- (struct sphor* sphor,
- double *out_estim,
- double *out_std);
+ (struct sphor* sphor);
#endif /* SPHOR_COMPUTE_MVREA_H */
