commit ae9d40be46e3d4d0e06d706261bc2db1da08dfe4
parent 722272532c4d654562730ce686b88e90441ea6f7
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Wed, 22 Oct 2025 16:43:30 +0200
Write estimators to output file
The introduction of accums and accum lists in sphor enables the
estimation of multiple random variables within a single algorithm. Each
random variable corresponds to an estimator that must be written to an
output file, either specified by the user or directed to the standard
output.
This commit implements this functionality by iterating over all accums
produced by the algorithm and, when necessary, combining certain accums
to construct others. This is required for random variables — such as the
total MVREA in the MVREA algorithm — that do not have directly
accessible accums during a realization, but are instead obtained as sums
of multiple accums. All resulting estimators are then written to the
output file.
In addition, this commit extends the internal accum API with several
utility functions and introduces the darray_accum_list structure and
programming interface, which facilitate the management of lists of
darray_accum objects (for instance when doing multi-thread
computations).
Diffstat:
5 files changed, 545 insertions(+), 48 deletions(-)
diff --git a/src/sphor.c b/src/sphor.c
@@ -102,6 +102,7 @@ sphor_create
sphor->allocator = allocator;
sphor->nthreads = args->nthreads;
sphor->samples = args->samples;
+ sphor->output_filename = args->output_filename;
if (NULL == args->logger) {
sphor->logger = LOGGER_DEFAULT;
} else {
diff --git a/src/sphor_accum.c b/src/sphor_accum.c
@@ -122,25 +122,142 @@ error:
}
extern LOCAL_SYM void
-accum_threads_aggregate
- (size_t nthreads,
- struct darray_accum* accums_threads,
+accum_compute_estim
+ (const struct accum* accum,
+ struct estim* estim)
+{
+ double sum = 0;
+ double sum2 = 0;
+ double samples = 0;
+ double mean = 0;
+ double var = 0;
+ double std = 0;
+
+ ASSERT(NULL != accum);
+ ASSERT(NULL != estim);
+
+ samples = (double)accum->n_realizations;
+ sum = accum->sum;
+ sum2 = accum->sum2;
+
+ mean = sum/samples;
+ var = (sum2/samples - mean * mean) / (samples - 1);
+ std = sqrt(var);
+
+ estim->mean = mean;
+ estim->var = var;
+ estim->std = std;
+}
+
+extern LOCAL_SYM void
+write_accum_estim
+ (const struct accum* accum,
+ FILE* stream)
+{
+ struct estim estim = ESTIM_NULL;
+
+ ASSERT(NULL != accum);
+ ASSERT(NULL != stream);
+
+ accum_compute_estim(accum, &estim);
+
+ if(0 != strcmp("\0", str_cget(&accum->component))){
+ fprintf
+ (stream, "%s:%s:%s:%lf:%lf\n",
+ str_cget(&accum->observable),
+ str_cget(&accum->sensor),
+ str_cget(&accum->component),
+ estim.mean, estim.std);
+ } else if (0 != strcmp("\0", str_cget(&accum->sensor))){
+ fprintf
+ (stream, "%s:%s:%lf:%lf\n",
+ str_cget(&accum->observable),
+ str_cget(&accum->sensor),
+ estim.mean, estim.std);
+ } else {
+ fprintf
+ (stream, "%s:%lf:%lf\n",
+ str_cget(&accum->observable),
+ estim.mean, estim.std);
+ }
+}
+
+extern LOCAL_SYM void
+accum_add
+ (const struct accum* op0,
+ const struct accum* op1,
+ struct accum* result)
+{
+ ASSERT(NULL != op0);
+ ASSERT(NULL != op1);
+ ASSERT(NULL != result);
+
+ result->sum = op0->sum + op1->sum;
+ result->sum2 = op0->sum2 + op1->sum2;
+ result->n_realizations = op0->n_realizations + op1->n_realizations;
+}
+
+extern LOCAL_SYM void
+accums_add
+ (const struct darray_accum* op0,
+ const struct darray_accum* op1,
+ struct darray_accum* result)
+{
+ size_t i = 0;
+ size_t naccums = 0;
+
+ ASSERT(NULL != op0);
+ ASSERT(NULL != op1);
+ ASSERT(NULL != result);
+
+ naccums = darray_accum_size_get(result);
+ ASSERT(darray_accum_size_get(op0) == darray_accum_size_get(op1));
+ ASSERT(darray_accum_size_get(op0) == naccums);
+
+ FOR_EACH(i, 0, naccums){
+ accum_add(darray_accum_cdata_get(op0) + i,
+ darray_accum_cdata_get(op1) + i,
+ darray_accum_data_get(result) + i);
+ }
+}
+
+extern LOCAL_SYM void
+sum_accums
+ (const struct darray_accum* accums,
+ struct accum* tot_accum)
+{
+ size_t j = 0;
+ size_t accum_count = 0;
+
+ ASSERT(NULL != accums);
+ ASSERT(NULL != tot_accum);
+
+ accum_count = darray_accum_size_get(accums);
+
+ FOR_EACH(j, 0, accum_count) {
+ struct accum* accum = NULL;
+ accum = (struct accum*)darray_accum_cdata_get(accums) + j;
+ accum_add(tot_accum, accum, tot_accum);
+ }
+}
+
+extern LOCAL_SYM void
+sum_accum_lists
+ (const struct darray_accum_list* accum_lists,
struct darray_accum* accums)
{
- size_t i = 0, j = 0;
+ size_t list_len = 0;
+ size_t i;
- ASSERT(NULL != accums_threads);
+ ASSERT(NULL != accum_lists);
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;
- }
+ list_len = darray_accum_list_size_get(accum_lists);
+
+ FOR_EACH(i, 0, list_len) {
+ struct darray_accum* accums_i = NULL;
+ accums_i = (struct darray_accum*)darray_accum_list_cdata_get(accum_lists)+i;
+ accums_add(accums, accums_i, accums);
}
}
@@ -157,8 +274,8 @@ update_accum
ASSERT(NULL != accum2id);
ASSERT(NULL != accums);
-
- paccum_id = htable_accum2id_find((struct htable_accum2id *)accum2id, accum_key);
+ paccum_id = htable_accum2id_find
+ ((struct htable_accum2id *)accum2id, accum_key);
ASSERT(NULL != paccum_id);
accum = &darray_accum_data_get(accums)[*paccum_id];
diff --git a/src/sphor_accum.h b/src/sphor_accum.h
@@ -41,6 +41,14 @@ struct accum {
#define ACCUM_NULL__ {0}
static const struct accum ACCUM_NULL = ACCUM_NULL__;
+struct estim {
+ double mean;
+ double var;
+ double std;
+};
+#define ESTIM_NULL__ {0}
+static const struct estim ESTIM_NULL = ESTIM_NULL__;
+
enum sphor_sensor_type {
SPHOR_SENSOR_VOLUME,
SPHOR_SENSOR_SURFACE,
@@ -74,6 +82,23 @@ accum_copy_and_release
(struct accum* dst,
struct accum* src);
+extern LOCAL_SYM void
+accum_compute_estim
+ (const struct accum* accum,
+ struct estim* estim);
+
+extern LOCAL_SYM void
+write_accum_estim
+ (const struct accum* accum,
+ FILE* stream);
+
+/* Add two accums and return it in result */
+extern LOCAL_SYM void
+accum_add
+ (const struct accum* op0,
+ const struct accum* op1,
+ struct accum* result);
+
/* Generate darray_accum data type and API */
#define DARRAY_NAME accum
#define DARRAY_DATA struct accum
@@ -83,6 +108,15 @@ accum_copy_and_release
#define DARRAY_FUNCTOR_COPY_AND_RELEASE accum_copy_and_release
#include <rsys/dynamic_array.h>
+/* Generate darray of darray_accum data type and API */
+#define DARRAY_NAME accum_list
+#define DARRAY_DATA struct darray_accum
+#define DARRAY_FUNCTOR_INIT darray_accum_init
+#define DARRAY_FUNCTOR_RELEASE darray_accum_release
+#define DARRAY_FUNCTOR_COPY darray_accum_copy
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_accum_copy_and_release
+#include <rsys/dynamic_array.h>
+
static INLINE char
accum_key_eq
(const struct accum_key* a,
@@ -142,10 +176,25 @@ register_accum
char* component,
struct darray_accum *accums);
+/* Add all the elements (sum) of an accum list and return it in result */
+extern LOCAL_SYM void
+sum_accums
+ (const struct darray_accum* accums,
+ struct accum* result);
+
+/* Add the elements of same index in two different accum lists and return it in
+ * result */
+extern LOCAL_SYM void
+accums_add
+ (const struct darray_accum* op0,
+ const struct darray_accum* op1,
+ struct darray_accum* result);
+
+/* Add all the elements (sum) of same index of a list of darray_accums and
+ * return it in result */
extern LOCAL_SYM void
-accum_threads_aggregate
- (size_t nthreads,
- struct darray_accum* accums_threads,
+sum_accum_lists
+ (const struct darray_accum_list* accum_lists,
struct darray_accum* accums);
extern LOCAL_SYM void
diff --git a/src/sphor_c.h b/src/sphor_c.h
@@ -69,6 +69,9 @@ struct sphor {
size_t samples; /* number of realisations used during the Monte Carlo
integrating algorithm */
+ /* Output stream (can be NULL: stdout) */
+ char* output_filename;
+
/* Log */
struct logger* logger;
int verbose;
diff --git a/src/sphor_compute_mvrea.c b/src/sphor_compute_mvrea.c
@@ -92,6 +92,7 @@ setup_MVREA_accums_surfaces
if (RES_OK != res){ goto error; }
}
}
+
exit:
return res;
error:
@@ -159,6 +160,7 @@ setup_MVREA_accums_volumes
}
}
}
+
exit:
return res;
error:
@@ -238,6 +240,7 @@ volume_compute_total_ka
*ka += concentration * sigma_a;
}
+
exit:
return res;
error:
@@ -271,6 +274,7 @@ sample_abs_free_path_exp
else {
*s = FLT_MAX;
}
+
exit:
return res;
error:
@@ -396,6 +400,341 @@ error:
goto exit;
}
+/*******************************************************************************
+ * Write volumes
+ ******************************************************************************/
+static res_T
+MVREA_aggregate_vol_accums
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ size_t i_volume,
+ struct darray_accum* vol_accums)
+{
+ char* volume_name;
+ size_t j_proprad = 0;
+ size_t prop_rad_count = 0;
+ size_t* paccum_id = NULL;
+ struct accum* accum = NULL;
+ struct accum vol_accum = ACCUM_NULL;
+ struct accum_key accum_key = ACCUM_KEY_NULL;
+ struct sphin_volume* volume = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+ ASSERT(NULL != vol_accums);
+
+ accum_init(sphor->allocator, &vol_accum);
+
+ res = sphin_config_get_volume(sphor->config, i_volume, &volume);
+ if (RES_OK != res){ goto error; }
+
+ 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; }
+
+ str_set(&vol_accum.observable, "MVREA");
+ str_set(&vol_accum.sensor, volume_name);
+
+ FOR_EACH(j_proprad, 0, prop_rad_count){
+ accum_key.sensor_type = SPHOR_SENSOR_VOLUME;
+ accum_key.sensor_id = i_volume;
+ accum_key.component_id = j_proprad;
+ paccum_id = htable_accum2id_find
+ ((struct htable_accum2id *)accum2id, &accum_key);
+ ASSERT(NULL != paccum_id);
+
+ accum = (struct accum*)&darray_accum_cdata_get(accums)[*paccum_id];
+ accum_add(&vol_accum, accum, &vol_accum);
+ }
+ res = darray_accum_push_back(vol_accums, &vol_accum);
+ if (RES_OK != res){ goto error; }
+
+exit:
+ accum_release(&vol_accum);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+MVREA_aggregate_vols_accums
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ struct darray_accum* vol_accums)
+{
+ size_t volume_count = 0;
+ size_t i_volume = 0;
+ struct sphin_config* config = sphor->config;
+ struct sphin_volume* volume = NULL;
+ struct sphin_sensor_volume* sensor = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+ ASSERT(NULL != vol_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 = MVREA_aggregate_vol_accums
+ (sphor, accum2id, accums, i_volume, vol_accums);
+ if (RES_OK != res){ goto error; }
+ }
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+MVREA_write_volumes_outputs
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ size_t samples,
+ FILE* stream)
+{
+ size_t i = 0;
+ struct accum* accum = NULL;
+ struct accum tot_accum = ACCUM_NULL;
+ struct darray_accum vols_accums;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+ ASSERT(NULL != stream);
+
+ accum_init(sphor->allocator, &tot_accum);
+ darray_accum_init(sphor->allocator, &vols_accums);
+ res = MVREA_aggregate_vols_accums(sphor, accum2id, accums, &vols_accums);
+ if (RES_OK != res){ goto error; }
+
+ sum_accums(&vols_accums, &tot_accum);
+
+ res = str_set(&tot_accum.observable, "MVREA");
+ if (RES_OK != res){ goto error; }
+ tot_accum.n_realizations = samples;
+
+ write_accum_estim(&tot_accum, stream);
+
+ accum = darray_accum_data_get(&vols_accums);
+ FOR_EACH(i, 0, darray_accum_size_get(&vols_accums)){
+ accum += i;
+ accum->n_realizations = samples;
+ write_accum_estim(accum, stream);
+ }
+
+exit:
+ accum_release(&tot_accum);
+ darray_accum_release(&vols_accums);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Write surfaces
+ ******************************************************************************/
+static res_T
+MVREA_aggregate_surf_accums
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ size_t i_surface,
+ struct darray_accum* surf_accums)
+{
+ char* surface_name;
+ size_t j_proprad = 0;
+ size_t* paccum_id = NULL;
+ struct accum* accum = NULL;
+ struct accum surf_accum = ACCUM_NULL;
+ struct accum_key accum_key = ACCUM_KEY_NULL;
+ struct sphin_surface* surface = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+ ASSERT(NULL != surf_accums);
+
+ accum_init(sphor->allocator, &surf_accum);
+
+ res = sphin_config_get_surface(sphor->config, i_surface, &surface);
+ if (RES_OK != res){ goto error; }
+
+ res = sphin_surface_get_name(surface, &surface_name);
+ if (RES_OK != res){ goto error; }
+
+ str_set(&surf_accum.observable, "LOSSES");
+ str_set(&surf_accum.sensor, surface_name);
+
+ accum_key.sensor_type = SPHOR_SENSOR_VOLUME;
+ accum_key.sensor_id = i_surface;
+ accum_key.component_id = j_proprad;
+ paccum_id = htable_accum2id_find
+ ((struct htable_accum2id *)accum2id, &accum_key);
+ ASSERT(NULL != paccum_id);
+
+ accum = (struct accum*)&darray_accum_cdata_get(accums)[*paccum_id];
+ accum_add(&surf_accum, accum, &surf_accum);
+
+ res = darray_accum_push_back(surf_accums, &surf_accum);
+ if (RES_OK != res){ goto error; }
+
+exit:
+ accum_release(&surf_accum);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+MVREA_aggregate_surfs_accums
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ struct darray_accum* surfs_accums)
+{
+ 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);
+ ASSERT(NULL != surfs_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 = MVREA_aggregate_surf_accums
+ (sphor, accum2id, accums, i_surface, surfs_accums);
+ if (RES_OK != res){ goto error; }
+ }
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+MVREA_write_surfaces_outputs
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ size_t samples,
+ FILE* stream)
+{
+ size_t i = 0;
+ struct accum* accum = NULL;
+ struct accum tot_accum = ACCUM_NULL;
+ struct darray_accum surf_accums;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+ ASSERT(NULL != stream);
+
+ accum_init(sphor->allocator, &tot_accum);
+ darray_accum_init(sphor->allocator, &surf_accums);
+ res = MVREA_aggregate_surfs_accums(sphor, accum2id, accums, &surf_accums);
+ if (RES_OK != res){ goto error; }
+
+ sum_accums(&surf_accums, &tot_accum);
+
+ res = str_set(&tot_accum.observable, "LOSSES");
+ if (RES_OK != res){ goto error; }
+
+ tot_accum.n_realizations = samples;
+
+ write_accum_estim(&tot_accum, stream);
+
+ accum = darray_accum_data_get(&surf_accums);
+ FOR_EACH(i, 0, darray_accum_size_get(&surf_accums)){
+ accum += i;
+ accum->n_realizations = samples;
+ write_accum_estim(accum, stream);
+ }
+
+exit:
+ accum_release(&tot_accum);
+ darray_accum_release(&surf_accums);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+MVREA_write_outputs
+ (struct sphor* sphor,
+ const struct htable_accum2id* accum2id,
+ const struct darray_accum* accums,
+ size_t nfailures)
+{
+ FILE* stream = NULL;
+ size_t samples = 0;
+ size_t i = 0;
+ struct accum* accum;
+ res_T res = RES_OK;
+
+ ASSERT(NULL != sphor);
+ ASSERT(NULL != accum2id);
+ ASSERT(NULL != accums);
+
+ if (NULL != sphor->output_filename){
+ stream = fopen(sphor->output_filename, "w");
+ } else { stream = stdout; }
+
+ samples = sphor->samples - nfailures;
+
+ res = MVREA_write_volumes_outputs(sphor, accum2id, accums, samples, stream);
+ if (RES_OK != res){ goto error; }
+
+ res = MVREA_write_surfaces_outputs(sphor, accum2id, accums, samples, stream);
+ if (RES_OK != res){ goto error; }
+
+ accum = (struct accum*)darray_accum_cdata_get(accums);
+ FOR_EACH(i, 0, darray_accum_size_get(accums)){
+ accum += i;
+ accum->n_realizations = samples;
+ write_accum_estim(accum, stream);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+
static res_T
compute_MVREA_realization
(struct sphor* sphor,
@@ -449,6 +788,7 @@ compute_MVREA_realization
* position, the photon is lost and it counts zero in the weight. */
if (S3D_HIT_NONE(&hit)) { break; }
+ /* TODO: Retrive prop rads of current volume */
/* Sample the free path length to absorption from an exponential
* distribution with rate parameter ka */
res = sample_abs_free_path_exp
@@ -528,6 +868,7 @@ compute_MVREA_realization
if (RES_OK != res){ goto error; }
d3_set_f3(pos, attrib.value);
}
+
exit:
return res;
error:
@@ -550,9 +891,7 @@ sphor_compute_MVREA
struct ssp_rng **rngs = 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;
+ struct darray_accum_list accums_threads;
samples = sphor->samples;
nthreads = sphor->nthreads;
@@ -570,14 +909,16 @@ sphor_compute_MVREA
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; }
+ darray_accum_list_init(sphor->allocator, &accums_threads);
+ res = darray_accum_list_resize(&accums_threads, nthreads);
+ if (RES_OK != res){ 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);
+ struct darray_accum* accum_thread;
+ accum_thread = darray_accum_list_data_get(&accums_threads) + i;
+ res = darray_accum_copy(accum_thread, &accums);
if (RES_OK != res){ goto error; }
}
@@ -602,13 +943,15 @@ sphor_compute_MVREA
for(i=0; i<samples; i++) {
const int ithread = omp_get_thread_num();
res_T res_local = RES_OK;
+ struct darray_accum* accum_thread;
+ accum_thread = darray_accum_list_data_get(&accums_threads) + ithread;
/* Ignore the rest of the for loop if there is a fatal (not failure) error
* in the previous realizations */
if (RES_OK != res) continue;
res_local = compute_MVREA_realization
- (sphor, rngs[ithread], &accum2id, &accums_threads[ithread]);
+ (sphor, rngs[ithread], &accum2id, accum_thread);
if (RES_OK != res_local) {
/*Protect res and nfailures from concurrent write accesses*/
#pragma omp critical
@@ -626,29 +969,12 @@ sphor_compute_MVREA
if (res != RES_OK) { goto error; }
/* 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);
+ sum_accum_lists(&accums_threads, &accums);
+
+ res = MVREA_write_outputs(sphor, &accum2id, &accums, nfailures);
+ if (res != RES_OK) { goto error; }
exit:
- /* 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) {
@@ -659,6 +985,7 @@ exit:
}
/* Free local data structures accums and accum2id */
darray_accum_release(&accums);
+ darray_accum_list_release(&accums_threads);
htable_accum2id_release(&accum2id);
return res;
