commit f3126f6015cf5b9b5a98d468eedca058ddc325bb
parent 6c95f8e3f25e03a8b2e606cbb183cba1269f6b0c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 23 Sep 2026 11:46:36 +0200
Merge remote-tracking branch 'origin/feature_sln_stat' into develop
Fix compilation issues due to API breaks on develop
Diffstat:
| M | .gitignore | | | 1 | + |
| M | Makefile | | | 20 | ++++++++++++++------ |
| A | doc/sln-stat.1 | | | 168 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/sln_stat.c | | | 486 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
4 files changed, 669 insertions(+), 6 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -10,6 +10,7 @@ mixture.txt
sln-build
sln-get
sln-slab
+sln-stat
tags
tags
test_*
diff --git a/Makefile b/Makefile
@@ -84,7 +84,7 @@ libsln.o: $(OBJ)
################################################################################
# Utils
################################################################################
-UTIL_SRC = src/sln_build.c src/sln_get.c src/sln_slab.c
+UTIL_SRC = src/sln_build.c src/sln_get.c src/sln_slab.c src/sln_stat.c
UTIL_OBJ = $(UTIL_SRC:.c=.o)
UTIL_DEP = $(UTIL_SRC:.c=.d)
@@ -103,14 +103,14 @@ LDFLAGS_SLAB = $(LDFLAGS_EXE) $(LIBS_SLAB)
utils: library $(UTIL_DEP) .config
@$(MAKE) -fMakefile \
$$(for i in $(UTIL_DEP); do printf -- '-f%s\n' "$${i}"; done) \
- sln-build sln-get sln-slab
+ sln-build sln-get sln-slab sln-stat
$(UTIL_DEP) $(UTIL_OBJ): config.mk sln-local.pc
-
src/sln_build.d src/sln_build.o: src/sln_build.c
src/sln_get.d src/sln_get.o: src/sln_get.c
src/sln_slab.d src/sln_slab.o: src/sln_slab.c
+src/sln_stat.d src/sln_stat.o: src/sln_stat.c
sln-build: config.mk sln-local.pc src/sln_build.o $(LIBNAME)
$(CC) $(CFLAGS_UTIL) -o $@ src/sln_build.o $(LDFLAGS_UTIL)
@@ -121,20 +121,23 @@ sln-get: config.mk sln-local.pc src/sln_get.o $(LIBNAME)
sln-slab: config.mk sln-local.pc src/sln_slab.o $(LIBNAME)
$(CC) $(CFLAGS_SLAB) -o $@ src/sln_slab.o $(LDFLAGS_SLAB)
+sln-stat: config.mk sln-local.pc src/sln_stat.o $(LIBNAME)
+ $(CC) $(CFLAGS_SLAB) -o $@ src/sln_stat.o $(LDFLAGS_SLAB)
+
src/sln_build.d src/sln_get.d:
@$(CC) $(CFLAGS_UTIL) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
-src/sln_slab.d:
+src/sln_slab.d src/sln_stat.d:
@$(CC) $(CFLAGS_SLAB) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/sln_build.o src/sln_get.o:
$(CC) $(CFLAGS_UTIL) -c $(@:.o=.c) -o $@
-src/sln_slab.o:
+src/sln_slab.o src/sln_stat.o:
$(CC) $(CFLAGS_SLAB) -c $(@:.o=.c) -o $@
clean_utils:
- rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build sln-get sln-slab
+ rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build sln-get sln-slab sln-stat
################################################################################
# Installation
@@ -172,11 +175,13 @@ install: library pkg utils
install 755 "$(DESTDIR)$(BINPREFIX)" sln-build; \
install 755 "$(DESTDIR)$(BINPREFIX)" sln-get; \
install 755 "$(DESTDIR)$(BINPREFIX)" sln-slab; \
+ install 755 "$(DESTDIR)$(BINPREFIX)" sln-stat; \
install 644 "$(DESTDIR)$(LIBPREFIX)/pkgconfig" sln.pc; \
install 644 "$(DESTDIR)$(INCPREFIX)/star" src/sln.h; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-build.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-get.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-slab.1; \
+ install 644 "$(DESTDIR)$(MANPREFIX)/man1" doc/sln-stat.1; \
install 644 "$(DESTDIR)$(MANPREFIX)/man5" doc/sln-mixture.5; \
install 644 "$(DESTDIR)$(PREFIX)/share/doc/star-line" COPYING README.md
@@ -185,12 +190,14 @@ uninstall:
rm -f "$(DESTDIR)$(BINPREFIX)/sln-build"
rm -f "$(DESTDIR)$(BINPREFIX)/sln-get"
rm -f "$(DESTDIR)$(BINPREFIX)/sln-slab"
+ rm -f "$(DESTDIR)$(BINPREFIX)/sln-stat"
rm -f "$(DESTDIR)$(LIBPREFIX)/pkgconfig/sln.pc"
rm -f "$(DESTDIR)$(BINPREFIX)/sln"
rm -f "$(DESTDIR)$(INCPREFIX)/star/sln.h"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-build.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-get.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-slab.1"
+ rm -f "$(DESTDIR)$(MANPREFIX)/man1/sln-stat.1"
rm -f "$(DESTDIR)$(MANPREFIX)/man5/sln-mixture.5"
rm -f "$(DESTDIR)$(PREFIX)/share/doc/star-line/COPYING"
rm -f "$(DESTDIR)$(PREFIX)/share/doc/star-line/README.md"
@@ -203,6 +210,7 @@ lint:
mandoc -Tlint -Wwarning doc/sln-build.1
mandoc -Tlint -Wwarning doc/sln-get.1
mandoc -Tlint -Wwarning doc/sln-slab.1
+ mandoc -Tlint -Wwarning doc/sln-stat.1
mandoc -Tlint doc/sln-mixture.5
################################################################################
diff --git a/doc/sln-stat.1 b/doc/sln-stat.1
@@ -0,0 +1,168 @@
+.\" Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+.\" Copyright (C) 2026 Université de Lorraine
+.\" Copyright (C) 2022 Centre National de la Recherche Scientifique
+.\" Copyright (C) 2022 Université Paul Sabatier
+.\"
+.\" This program is free software: you can redistribute it and/or modify
+.\" it under the terms of the GNU General Public License as published by
+.\" the Free Software Foundation, either version 3 of the License, or
+.\" (at your option) any later version.
+.\"
+.\" This program is distributed in the hope that it will be useful,
+.\" but WITHOUT ANY WARRANTY; without even the implied warranty of
+.\" MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+.\" GNU General Public License for more details.
+.\"
+.\" You should have received a copy of the GNU General Public License
+.\" along with this program. If not, see <http://www.gnu.org/licenses/>.
+.Dd May 5, 2026
+.Dt SLN-STAT 1
+.Os
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh NAME
+.Nm sln-stat
+.Nd computations of basic statistics over k spectrum
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh SYNOPSIS
+.Nm
+.Op Fl dhsv
+.Op Fl n Ar nrealisations
+.Op Fl t Ar threads
+.Fl S Ar nu_min , Ns Ar nu_max
+.Fl a Ar accel_struct
+.Fl m Ar molparams
+.Fl l Ar lines
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh DESCRIPTION
+.Nm
+calculates the mean and mean of squares of k absorption values using a Monte
+Carlo algorithm that samples the spectral lines that make up the gas mixture.
+These computations are accelerated by sampling the lines based on the magnitude
+of their contribution to the mixture’s spectrum, so that few Monte Carlo runs
+are required to estimate the spectrum statistics with a high degree of
+confidence.
+The core of the proposal rests on this sampling strategy, made possible
+by constructing an acceleration structure from the set of lines in the
+mixture.
+A structure built using the
+.Xr sln-build 1
+utility and provided as input to the program.
+.Pp
+The output of
+.Nm
+displays the estimated mean and mean of squares, their standard
+deviation, and the number of Monte Carlo realisations rejected due to
+issues encountered during the computation, such as numerical
+uncertainty.
+Each estimate is displayed on a line formatted as follows:
+.Bd -literal -offset Ds
+"%-16s: %e +/- %e; %lu\en", name, estimate, std_err, rejects_count
+.Ed
+.Pp
+The options are as follows:
+.Bl -tag -width Ds
+.\""""""""""""""""""""""""""""""""""
+.It Fl a Ar accel_struct
+An acceleration structure corresponding to the input
+.Ar lines ,
+used to accelerate their sampling based on their importance.
+This structure is generated by the
+.Xr sln-build 1
+tool.
+.\""""""""""""""""""""""""""""""""""
+.It Fl d
+Disables verification of the correspondence between the lines provided
+by the
+.Fl l
+option and those used to construct
+the acceleration structure defined by the
+.Fl a
+option.
+.Pp
+Warning!
+It is always recommended to verify that the data is correct, even though
+this verification can take a significant amount of time when there are a
+large number of lines.
+Anyway, a user who is
+.Em certain
+of the data’s consistency may nevertheless use this option
+.Pq at their own risk
+to disable this verification and thus speed up the execution.
+.\""""""""""""""""""""""""""""""""""
+.It Fl h
+Display short help and exit.
+.It Fl l Ar lines
+List of lines from which the tree was built.
+This list is in binary format as generated by the
+.Xr shtr 1
+binary, or in plain text HITRAN format, depending on whether the
+.Fl s
+option is set or not, respectively.
+.\""""""""""""""""""""""""""""""""""
+.It Fl m Ar molparams
+Isotopologue metadata in HITRAN format.
+.\""""""""""""""""""""""""""""""""""
+.It Fl n Ar nrealisations
+Number of Monte Carlo realisations.
+By default the number of realisations is 10000.
+.\""""""""""""""""""""""""""""""""""
+.It Fl S Ar nu_min , Ns Ar nu_max
+The spectral range, in cm^-1, over which the computations are performed.
+The default spectral range is from 0 to infinity.
+.\""""""""""""""""""""""""""""""""""
+.It Fl s
+Specifies that input lines are formatted according to the binary format
+as written by the
+.Xr shtr 1
+utility, and not according to the HITRAN format.
+This format is more compact, allowing for faster loading of line data.
+.\""""""""""""""""""""""""""""""""""
+.It Fl t Ar threads
+Advice on the number of threads to use.
+By default,
+.Nm
+uses as many threads as processor cores.
+.\""""""""""""""""""""""""""""""""""
+.It Fl v
+Make
+.Nm
+verbose.
+Multiple
+.Fl v
+options increase the verbosity.
+The maximum is 3.
+.El
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh EXIT STATUS
+.Ex -std
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh EXAMPLES
+Estimate the mean k and mean square k between 100 and 2500 cm^-1
+for a gaz mixture made of H2O, CO2 and CO molecules.
+The thermodynamic properties of the mixture, such as its pressure,
+temperature and molecular concentrations, correspond to those used to
+construct the acceleration structures with sln-build, provided as input
+arguments
+.Pq option Fl a .
+The isotopic metadata
+.Pq option Fl m
+and the list of lines
+.Pq option Fl l
+partitioned by the acceleration structure, complete the list of input
+data.
+The latter is encoded in the format generated by the
+.Xr shtr 1
+tool
+.Pq option Fl s .
+The isotopes are in HITRAN format.
+Finally, make the program as verbose as possible
+.Pq options Fl vvv .
+.Bd -literal -offset Ds
+sln-stat -S 100,2500 -a tree_H2O_CO2_CO_1atm_600K.sln \e
+ -m molparam.txt -sl H2O_CO2_CO_100-2500cm-1.shtr -vvv
+.Ed
+.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+.Sh SEE ALSO
+.Xr shtr 1 ,
+.Xr sln-build 1 ,
+.Xr sln-slab 1
diff --git a/src/sln_stat.c b/src/sln_stat.c
@@ -0,0 +1,486 @@
+/* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2026 Université de Lorraine
+ * Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* getopt */
+
+#include "sln.h"
+
+#include <star/shtr.h>
+#include <star/sbb.h>
+#include <star/ssp.h>
+
+#include <rsys/cstr.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/str.h>
+
+#include <omp.h>
+
+#include <unistd.h> /* getopt */
+
+enum estimate {
+ MEAN,
+ SQMEAN,
+ ESTIMATE_COUNT__
+};
+
+#define WAVENUMBER_TO_WAVELENGTH(Nu/* [cm^-1] */) (1.e-2/(Nu))/*[m]*/
+
+struct args {
+ const char* tree; /* Acceleration structure */
+ const char* molparams;
+ const char* lines;
+
+ double spectral_range[2]; /* [cm^-1]^2 */
+
+ unsigned long nrealisations; /* Number of Monte Carlo realisations */
+
+ /* Miscellaneous */
+ unsigned nthreads_hint; /* Hint on the number of threads to use */
+ int disable_line_hash_check;
+ int lines_in_shtr_format;
+ int verbose;
+ int quit;
+};
+#define ARGS_DEFAULT__ {NULL,NULL,NULL,{0,DBL_MAX},10000,UINT_MAX,0,0,0,0}
+static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
+
+struct cmd {
+ struct args args;
+
+ struct sln_tree* tree;
+ unsigned nthreads;
+};
+#define CMD_NULL__ {0}
+static const struct cmd CMD_NULL = CMD_NULL__;
+
+struct accum {
+ double sum;
+ double sum2;
+ size_t count;
+};
+#define ACCUM_NULL__ {0}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+usage(FILE* stream)
+{
+ fprintf(stream,
+"usage: sln-stat [-dhsv] [-n nrealisations] [-t threads]\n"
+" -S nu_min,nu_max -a accel_struct -m molparams -l lines\n");
+}
+
+static res_T
+parse_spectral_range(const char* str, double spectral_range[2])
+{
+ size_t len = 0;
+ res_T res = RES_OK;
+ ASSERT(str && spectral_range);
+
+ res = cstr_to_list_double(str, ',', spectral_range, &len, 2);
+ if(res == RES_OK && len < 2) res = RES_BAD_ARG;
+
+ return res;
+}
+
+static res_T
+args_init(struct args* args, int argc, char** argv)
+{
+ int opt = 0;
+ res_T res = RES_OK;
+
+ ASSERT(args);
+
+ *args = ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "a:dhl:m:n:S:st:v")) != -1) {
+ switch(opt) {
+ case 'a': args->tree = optarg; break;
+ case 'd': args->disable_line_hash_check = 1; break;
+ case 'h':
+ usage(stdout);
+ args->quit = 1;
+ goto exit;
+ case 'l': args->lines = optarg; break;
+ case 'm': args->molparams = optarg; break;
+ case 'n': res = cstr_to_ulong(optarg, &args->nrealisations); break;
+ case 'S': res = parse_spectral_range(optarg, args->spectral_range); break;
+ case 's': args->lines_in_shtr_format = 1; break;
+ case 't':
+ res = cstr_to_uint(optarg, &args->nthreads_hint);
+ if(res == RES_OK && args->nthreads_hint == 0) res = RES_BAD_ARG;
+ break;
+ case 'v': args->verbose += (args->verbose < 3); break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ #define MANDATORY(Cond, Name, Opt) { \
+ if(!(Cond)) { \
+ fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \
+ res = RES_BAD_ARG; \
+ goto error; \
+ } \
+ } (void)0
+ MANDATORY(args->molparams, "molparams", 'm');
+ MANDATORY(args->lines, "line list", 'l');
+ MANDATORY(args->tree, "acceleration structure", 'a');
+ #undef MANDATORY
+
+exit:
+ return res;
+error:
+ usage(stderr);
+ goto exit;
+}
+
+static res_T
+load_lines
+ (struct shtr* shtr,
+ const struct args* args,
+ struct shtr_line_list** out_lines)
+{
+ struct shtr_line_list* lines = NULL;
+ res_T res = RES_OK;
+ ASSERT(shtr && args && out_lines);
+
+ if(args->lines_in_shtr_format) {
+ struct shtr_line_list_read_args read_args = SHTR_LINE_LIST_READ_ARGS_NULL;
+
+ /* Loads lines from data serialized by the Star-HITRAN library */
+ read_args.filename = args->lines;
+ res = shtr_line_list_read(shtr, &read_args, &lines);
+ if(res != RES_OK) goto error;
+
+ } else {
+ struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
+
+ /* Loads lines from a file in HITRAN format */
+ load_args.filename = args->lines;
+ res = shtr_line_list_load(shtr, &load_args, &lines);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_lines = lines;
+ return res;
+error:
+ if(lines) { SHTR(line_list_ref_put(lines)); lines = NULL; }
+ goto exit;
+}
+
+static void
+delete_per_thread_rngs(const struct cmd* cmd, struct ssp_rng* rngs[])
+{
+ unsigned i = 0;
+ ASSERT(cmd && rngs);
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ mem_rm(rngs);
+}
+
+static res_T
+create_per_thread_rngs(const struct cmd* cmd, struct ssp_rng** out_rngs[])
+{
+ struct ssp_rng_proxy* proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t i = 0;
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ rngs = mem_calloc(cmd->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+
+ res = ssp_rng_proxy_create(NULL, SSP_RNG_THREEFRY, cmd->nthreads, &proxy);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ res = ssp_rng_proxy_create_rng(proxy, i, &rngs[i]);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_rngs = rngs;
+ if(proxy) SSP(rng_proxy_ref_put(proxy));
+ return res;
+error:
+ if(cmd->args.verbose >= 1) {
+ fprintf(stderr,
+ "Error creating the list of per thread RNG -- %s\n",
+ res_to_cstr(res));
+ }
+ if(rngs) delete_per_thread_rngs(cmd, rngs);
+ rngs = NULL;
+ goto exit;
+}
+
+static void
+cmd_release(struct cmd* cmd)
+{
+ ASSERT(cmd);
+ if(cmd->tree) SLN(tree_ref_put(cmd->tree));
+}
+
+static res_T
+cmd_init(struct cmd* cmd, const struct args* args)
+{
+ /* Star Line */
+ struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
+ struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL;
+ struct sln_device* sln = NULL;
+
+ /* Star HITRAN */
+ struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
+ struct shtr* shtr = NULL;
+ struct shtr_isotope_metadata* molparams = NULL;
+ struct shtr_line_list* lines = NULL;
+
+ /* Miscellaneous */
+ unsigned nthreads_max = 0;
+ res_T res = RES_OK;
+
+ ASSERT(cmd && args);
+
+ *cmd = CMD_NULL;
+
+ shtr_args.verbose = args->verbose;
+ res = shtr_create(&shtr_args, &shtr);
+ if(res != RES_OK) goto error;
+
+ res = shtr_isotope_metadata_load(shtr, args->molparams, &molparams);
+ if(res != RES_OK) goto error;
+
+ res = load_lines(shtr, args, &lines);
+ if(res != RES_OK) goto error;
+
+ sln_args.verbose = args->verbose;
+ res = sln_device_create(&sln_args, &sln);
+ if(res != RES_OK) goto error;
+
+ tree_args.metadata = molparams;
+ tree_args.lines = lines;
+ tree_args.filename = args->tree;
+ tree_args.disable_line_hash_check = args->disable_line_hash_check;
+ res = sln_tree_read(sln, &tree_args, &cmd->tree);
+ if(res != RES_OK) goto error;
+
+ nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs());
+ cmd->args = *args;
+ cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max);
+
+exit:
+ if(sln) SLN(device_ref_put(sln));
+ if(shtr) SHTR(ref_put(shtr));
+ if(molparams) SHTR(isotope_metadata_ref_put(molparams));
+ if(lines) SHTR(line_list_ref_put(lines));
+ return res;
+error:
+ cmd_release(cmd);
+ *cmd = CMD_NULL;
+ goto exit;
+}
+
+static INLINE const char*
+estimate_cstr(const enum estimate estimate)
+{
+ const char* cstr = NULL;
+ switch(estimate) {
+ case MEAN: cstr="mean"; break;
+ case SQMEAN: cstr="mean_of_squares"; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return cstr;
+}
+
+static res_T
+realisation
+ (const struct cmd* cmd,
+ struct ssp_rng* rng,
+ double out_weights[ESTIMATE_COUNT__])
+{
+ /* Acceleration structure */
+ struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL;
+ const struct sln_node* root = NULL;
+
+ /* Variables to sample k */
+ const struct sln_node* leaf1 = NULL;
+ const struct sln_node* leaf2 = NULL;
+ double leaf_proba1 = 0; /* Probability of sampling a line */
+ double leaf_proba2 = 0; /* Probability of sampling a line */
+ double leaf_ka1 = 0; /* Value of a line */
+ double leaf_ka2 = 0; /* Value of a line */
+
+ /* Miscellaneous */
+ double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */
+ double nu = 0; /* Sampled wavenumber [cm^-1] */
+ int i = 0;
+ res_T res = RES_OK;
+
+ ASSERT(cmd && rng && out_weights); /* Check pre-conditions */
+
+ /* Uniformly sample the spectral dimension */
+ nu = ssp_rng_uniform_double
+ (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]);
+
+ SLN(tree_get_desc(cmd->tree, &tree_desc));
+
+ /* Store the root node of the tree */
+ root = sln_tree_get_root(cmd->tree);
+
+ /* Importance sampling of a line and evaluation of the line contribution */
+ leaf1 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba1);
+ if(!leaf1) { res = RES_BAD_ARG; goto error; }
+ leaf_ka1 = sln_node_eval(cmd->tree, leaf1, NULL, nu);
+
+ /* Importance sampling of a line and evaluation of the line contribution */
+ leaf2 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba2);
+ if(!leaf2) { res = RES_BAD_ARG; goto error; }
+ leaf_ka2 = sln_node_eval(cmd->tree, leaf2, NULL, nu);
+
+ w[MEAN] = leaf_ka1 / leaf_proba1;
+ w[SQMEAN] = (leaf_ka1/leaf_proba1)*(leaf_ka2/leaf_proba2);
+
+exit:
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i];
+ return res;
+error:
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN;
+ goto exit;
+}
+
+static res_T
+cmd_run(const struct cmd* cmd)
+{
+ /* Random Number Generator */
+ struct ssp_rng** rngs = NULL;
+
+ /* Monte Carlo */
+ struct accum accum[ESTIMATE_COUNT__] = {0};
+ int64_t i = 0; /* Index of the realisation */
+ size_t nrejects = 0; /* Number of rejected realisations */
+
+ /* Progress */
+ size_t nrealisations = 0;
+ size_t realisation_done = 0;
+ int progress = 0;
+ int progress_pcent = 10;
+
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ res = create_per_thread_rngs(cmd, &rngs);
+ if(res != RES_OK) goto error;
+
+ #define PROGRESS_MSG "Solving: %3d%%\n"
+ if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress);
+
+ nrealisations = cmd->args.nrealisations;
+
+ omp_set_num_threads((int)cmd->nthreads);
+
+ #pragma omp parallel for schedule(static)
+ for(i = 0; i < (int64_t)nrealisations; ++i) {
+ double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */
+ const int ithread = omp_get_thread_num();
+ int pcent = 0;
+ res_T res_realisation = RES_OK;
+
+ res_realisation = realisation(cmd, rngs[ithread], w);
+
+ #pragma omp critical
+ {
+ /* Update the Monte Carlo accumulator */
+ if(res_realisation == RES_OK) {
+ int iestim = 0;
+ FOR_EACH(iestim, 0, ESTIMATE_COUNT__) {
+ accum[iestim].sum += w[iestim];
+ accum[iestim].sum2 += w[iestim]*w[iestim];
+ accum[iestim].count += 1;
+ }
+ }
+
+ if(cmd->args.verbose >= 3) {
+ /* Update progress */
+ realisation_done += 1;
+ pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5);
+ if(pcent/progress_pcent > progress/progress_pcent) {
+ progress = pcent;
+ fprintf(stderr, PROGRESS_MSG, progress);
+ }
+ }
+ }
+ }
+
+ #undef PROGRESS_MSG
+
+ nrejects = nrealisations - accum[0].count;
+
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) {
+ const double E = accum[i].sum / (double)accum[i].count;
+ const double V = accum[i].sum2 / (double)accum[i].count - E*E;
+ const double SE = sqrt(V/(double)accum[i].count);
+
+ /* Assume that the number of realisations is the same for all estimates */
+ ASSERT(accum[i].count == accum[0].count);
+
+ printf("%-16s: %e +/- %e; %lu\n",
+ estimate_cstr(i), E, SE, (unsigned long)nrejects);
+ }
+
+exit:
+ delete_per_thread_rngs(cmd, rngs);
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct args args = ARGS_DEFAULT;
+ struct cmd cmd = CMD_NULL;
+ int err = 0;
+ res_T res = RES_OK;
+
+ if((res = args_init(&args, argc, argv)) != RES_OK) goto error;
+ if(args.quit) goto exit;
+
+ if((res = cmd_init(&cmd, &args)) != RES_OK) goto error;
+ if((res = cmd_run(&cmd)) != RES_OK) goto error;
+
+exit:
+ cmd_release(&cmd);
+ CHK(mem_allocated_size() == 0);
+ return err;
+error:
+ err = 1;
+ goto exit;
+}
