commit 8a1a290654be84a475518257e6473aea5bde8005
parent 305f4f4329e80e77b3906141203e9241be7188c0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 12 May 2026 12:39:49 +0200
Add the sln_node_sample_leaf function
It implements the importance sampling procedure proposed by
Nyffenegger-Péré et al. [1], which was previously implemented in the
sln-slab tool. This means that this procedure is sufficiently
well-defined to be provided directly as an API function. Of course, the
caller can still implement their own original sampling function by
explicitly traversing the tree, as sln-slab did before the library
directly provided the sampling function it used.
[1] Y. Nyffenegger-Péré et al., “Spectrally refined unbiased monte
carlo estimate of the earth's global radiative cooling”, Proceedings
of the National Academy of Sciences, 5, 121, e2315492121, 2024.
Diffstat:
6 files changed, 149 insertions(+), 127 deletions(-)
diff --git a/Makefile b/Makefile
@@ -143,8 +143,9 @@ pkg:
sed -e 's#@PREFIX@#$(PREFIX)#g' \
-e 's#@VERSION@#$(VERSION)#g' \
-e 's#@LBLU_VERSION@#$(LBLU_VERSION)#g' \
- -e 's#@SHTR_VERSION@#$(SHTR_VERSION)#g' \
-e 's#@RSYS_VERSION@#$(RSYS_VERSION)#g' \
+ -e 's#@SHTR_VERSION@#$(SHTR_VERSION)#g' \
+ -e 's#@SSP_VERSION@#$(SSP_VERSION)#g' \
sln.pc.in > sln.pc
sln-local.pc: sln.pc.in
@@ -153,8 +154,9 @@ sln-local.pc: sln.pc.in
-e 's#^libdir=.*#libdir=./#' \
-e 's#@VERSION@#$(VERSION)#g' \
-e 's#@LBLU_VERSION@#$(LBLU_VERSION)#g' \
- -e 's#@SHTR_VERSION@#$(SHTR_VERSION)#g' \
-e 's#@RSYS_VERSION@#$(RSYS_VERSION)#g' \
+ -e 's#@SHTR_VERSION@#$(SHTR_VERSION)#g' \
+ -e 's#@SSP_VERSION@#$(SSP_VERSION)#g' \
sln.pc.in > $@
install: library pkg utils
@@ -210,13 +212,14 @@ TEST_SRC =\
src/test_sln_device.c\
src/test_sln_mesh.c\
src/test_sln_mixture.c\
- src/test_sln_tree.c
+ src/test_sln_tree.c\
+ src/test_sln_tree_sample.c
TEST_OBJ = $(TEST_SRC:.c=.o)
TEST_DEP = $(TEST_SRC:.c=.d)
TEST_TGT = $(TEST_SRC:.c=.t)
-INCS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local shtr)
-LIBS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr)
+INCS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local shtr star-sp)
+LIBS_TEST = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr star-sp)
CFLAGS_TEST = -std=c89 $(CFLAGS_EXE) $(INCS_TEST)
LDFLAGS_TEST = $(LDFLAGS_EXE) $(LIBS_TEST)
@@ -243,6 +246,7 @@ $(TEST_OBJ): config.mk sln-local.pc
test_sln_device \
test_sln_mesh \
test_sln_mixture \
+test_sln_tree_sample \
: config.mk sln-local.pc $(LIBNAME)
$(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST)
diff --git a/config.mk b/config.mk
@@ -31,13 +31,13 @@ PCFLAGS = $(PCFLAGS_$(LIB_TYPE))
LBLU_VERSION = 0.0.1
RSYS_VERSION = 0.14
SHTR_VERSION = 0
+SSP_VERSION = 0.15
# Dependencies for the sln-slab util
-SSP_VERSION = 0.15
SBB_VERSION = 0.0
-INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags lblu rsys shtr) -fopenmp
-LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs lblu rsys shtr) -fopenmp -lm
+INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags lblu rsys shtr star-sp) -fopenmp
+LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs lblu rsys shtr star-sp) -fopenmp -lm
################################################################################
# Compilation options
diff --git a/sln.pc.in b/sln.pc.in
@@ -3,7 +3,10 @@ includedir=${prefix}/include
libdir=${prefix}/lib
Requires: rsys >= @RSYS_VERSION@
-Requires.private: lblu >= @LBLU_VERSION@, shtr >= @SHTR_VERSION@
+Requires.private:\
+ lblu >= @LBLU_VERSION@,\
+ shtr >= @SHTR_VERSION@, \
+ star-sp >= @SSP_VERSION@
Name: Star-Line
Description: Star Line library
Version: @VERSION@
diff --git a/src/sln.h b/src/sln.h
@@ -254,6 +254,9 @@ struct sln_line {
#define SLN_LINE_NULL__ {0,0,0,0,SHTR_MOLECULE_ID_NULL}
static const struct sln_line SLN_LINE_NULL = SLN_LINE_NULL__;
+/* External data structure */
+struct ssp_rng;
+
/* Forward declarations of opaque data structures */
struct sln_device;
struct sln_mixture;
@@ -392,6 +395,16 @@ sln_node_get_mesh
const struct sln_node* node,
struct sln_mesh* mesh);
+/* Sample a leaf based on its importance, that is, its contribution to the
+ * node's absorption spectrum at the given wave number */
+SLN_API const struct sln_node* /* NULL <=> an error occurs */
+sln_node_sample_leaf
+ (const struct sln_tree* tree,
+ const struct sln_node* node,
+ const double nu, /* [cm^-1] */
+ struct ssp_rng* rng,
+ double* proba); /* May be NULL */
+
/*******************************************************************************
* Miscellaneous
******************************************************************************/
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -311,118 +311,6 @@ error:
goto exit;
}
-static res_T
-build_node_cumulative
- (const struct cmd* cmd,
- const struct sln_node* node,
- const double nu, /* [cm^-1] */
- const char* path, /* String representing the path to the node */
- double probabilities[SLN_TREE_ARITY_MAX],
- double cumulative[SLN_TREE_ARITY_MAX])
-{
- struct sln_mesh mesh = SLN_MESH_NULL;
- double ka = 0;
- double cumul = 0;
- unsigned i=0, n=0;
- res_T res = RES_OK;
- ASSERT(cmd && node && path && cumulative);
-
- n = sln_node_get_child_count(cmd->tree, node);
- ASSERT(n <= SLN_TREE_ARITY_MAX);
-
- FOR_EACH(i, 0, n) {
- const struct sln_node* child = sln_node_get_child(cmd->tree, node, i);
-
- SLN(node_get_mesh(cmd->tree, child, &mesh));
-
- ka = sln_mesh_eval(&mesh, nu);
-
- cumul += ka;
- cumulative[i] = cumul;
- probabilities[i] = ka;
- }
-
- /* Check the criterion of transition importance sampling, i.e. the value of
- * the parent node is greater than or equal to the sum of the values of its
- * children */
- SLN(node_get_mesh(cmd->tree, node, &mesh));
- ka = sln_mesh_eval(&mesh, nu);
- if(ka < cumul) {
- if(cmd->args.verbose >= 1) {
- fprintf(stderr,
- "error: ka < ka_{0} + ka_{1} + ... ka_{N-1}; "
- "%e < %e; nu=%-21.20g cm^-1; node path: %s\n",
- ka, cumul, nu, path);
- }
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Complete the probability calculation and normalize the cumulative */
- ASSERT(cumul != 0);
- FOR_EACH(i, 0, n) probabilities[i] /= cumul;
- FOR_EACH(i, 0, n-1) cumulative[i] /= cumul;
- cumulative[n-1] = 1; /* Handle numerical uncertainty */
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static const struct sln_node* /* NULL <=> an error occurs */
-sample_leaf
- (const struct cmd* cmd,
- struct ssp_rng* rng,
- const double nu/*[cm^-1]*/,
- double* out_proba)
-{
- double cumulative[SLN_TREE_ARITY_MAX];
- double probabilities[SLN_TREE_ARITY_MAX];
- struct str path;
- const struct sln_node* node = NULL;
- double proba = 1; /* Proba to sample the line */
- size_t depth = 0;
- res_T res = RES_OK;
- ASSERT(cmd && rng);
-
- str_init(NULL, &path);
- node = sln_tree_get_root(cmd->tree);
-
- for(depth=0; !sln_node_is_leaf(node); ++depth) {
- double r = 0;
- unsigned ichild = 0;
-
- res = build_node_cumulative
- (cmd, node, nu, str_cget(&path), probabilities, cumulative);
- if(res != RES_OK) goto error;
-
- /* Sample a child with respect to its importance */
- r = ssp_rng_canonical(rng);
- FOR_EACH(ichild, 0, SLN_TREE_ARITY_MAX) {
- if(r < cumulative[ichild]) {
- proba *= probabilities[ichild];
- node = sln_node_get_child(cmd->tree, node, ichild);
- break;
- }
-
- /* Update the path string */
- res = str_append_printf(&path, " -c%d", ichild);
- if(res != RES_OK) goto error;
- }
- ASSERT(ichild < SLN_TREE_ARITY_MAX); /* A node should have been sampled */
- }
-
-exit:
- str_release(&path);
- *out_proba = proba;
- return node;
-error:
- proba = NaN;
- node = NULL;
- goto exit;
-}
-
/* Check that the probability is valid */
static INLINE res_T
check_proba(const struct cmd* cmd, const double proba)
@@ -467,7 +355,7 @@ realisation
{
/* Acceleration structure */
struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL;
- const struct sln_node* node = NULL;
+ const struct sln_node* root = NULL;
struct sln_mesh mesh = SLN_MESH_NULL;
/* Null collisions */
@@ -497,11 +385,12 @@ realisation
SLN(tree_get_desc(cmd->tree, &tree_desc));
/* Retrieve the ka_max of the spectrum at the sampled nu */
- node = sln_tree_get_root(cmd->tree);
- SLN(node_get_mesh(cmd->tree, node, &mesh));
+ root = sln_tree_get_root(cmd->tree);
+ SLN(node_get_mesh(cmd->tree, root, &mesh));
ka_max = sln_mesh_eval(&mesh, nu);
for(ncollisions=0; ; ++ncollisions) {
+ const struct sln_node* leaf = NULL;
double dst = 0; /* Sampled distance */
double proba_abs = 0; /* Probability of absorption */
double leaf_proba = 0; /* Probability of sampling the line */
@@ -517,12 +406,12 @@ realisation
}
/* Importance sampling of a line */
- node = sample_leaf(cmd, rng, nu, &leaf_proba);
- if(!node) { res = RES_BAD_ARG; goto error; }
+ leaf = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba);
+ if(!leaf) { res = RES_BAD_ARG; goto error; }
/* Evaluate the value of the line and compute the probability of being
* absorbed by it */
- leaf_ka = sln_node_eval(cmd->tree, node, nu);
+ leaf_ka = sln_node_eval(cmd->tree, leaf, nu);
proba_abs = leaf_ka / (leaf_proba*ka_max);
if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error;
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -22,6 +22,7 @@
#include "sln_tree_c.h"
#include <star/shtr.h>
+#include <star/ssp.h>
#include <rsys/algorithm.h>
#include <rsys/cstr.h>
@@ -382,6 +383,67 @@ cmp_nu_vtx(const void* key, const void* item)
return 0;
}
+static res_T
+build_node_cumulative
+ (const struct sln_tree* tree,
+ const struct sln_node* node,
+ const double nu, /* [cm^-1] */
+ double proba[SLN_TREE_ARITY_MAX],
+ double cumul[SLN_TREE_ARITY_MAX])
+{
+ struct sln_mesh mesh = SLN_MESH_NULL;
+ double ka = 0;
+ double sum = 0;
+ unsigned i=0, n=0;
+ res_T res = RES_OK;
+ ASSERT(tree && node && proba && cumul);
+
+ n = sln_node_get_child_count(tree, node);
+ ASSERT(n <= SLN_TREE_ARITY_MAX);
+
+ FOR_EACH(i, 0, n) {
+ const struct sln_node* child = sln_node_get_child(tree, node, i);
+
+ SLN(node_get_mesh(tree, child, &mesh));
+ ka = sln_mesh_eval(&mesh, nu);
+
+ sum += ka;
+ cumul[i] = sum;
+ proba[i] = ka;
+ }
+
+ /* No lines could be sampled because none of them affect the absorption at the
+ * given wave number */
+ if(sum == 0) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check the criterion of transition importance sampling, i.e. the value of
+ * the parent node must be greater than or equal to the sum of the values of
+ * its children */
+ SLN(node_get_mesh(tree, node, &mesh));
+ ka = sln_mesh_eval(&mesh, nu);
+ if(ka < sum) {
+ ERROR(tree->sln,
+ "ka < ka_{0} + ka_{1} + ... + ka_{N-1}; %e < %e; nu=%-21.20g cm^-1\n",
+ ka, sum, nu);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Complete the probability calculation and normalize the cumulative */
+ ASSERT(sum != 0);
+ FOR_EACH(i, 0, n) proba[i] /= sum;
+ FOR_EACH(i, 0, n-1) cumul[i] /= sum;
+ cumul[n-1] = 1; /* Handle numerical uncertainty */
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
static void
release_tree(ref_T* ref)
{
@@ -749,6 +811,57 @@ sln_node_get_desc
return RES_OK;
}
+const struct sln_node*
+sln_node_sample_leaf
+ (const struct sln_tree* tree,
+ const struct sln_node* root,
+ const double nu, /*[cm^-1]*/
+ struct ssp_rng* rng,
+ double* out_leaf_proba) /* May be NULL */
+{
+ /* Temporary buffers used to store the cumulative of child nodes and their
+ * probability of being sampled based on their importance */
+ double cumul[SLN_TREE_ARITY_MAX];
+ double proba[SLN_TREE_ARITY_MAX];
+
+ const struct sln_node* node = NULL;
+ double leaf_proba = 1;
+ int depth = 0;
+ res_T res = RES_OK;
+
+ if(!tree || !root || !rng) { res = RES_BAD_ARG; goto error; }
+
+ for(depth=0, node=root; !sln_node_is_leaf(node); ++depth) {
+ double r = 0; /* Random number */
+ unsigned ichild = 0;
+
+ res = build_node_cumulative(tree, node, nu, proba, cumul);
+ if(res != RES_OK) goto error;
+
+ /* Sample a child node based on its importance. Use a simple linear search,
+ * since the tree's arity is small enough that a binary search is not
+ * necessary. FIXME if performance measurements show that this linear search
+ * incurs a significant cost */
+ r = ssp_rng_canonical(rng);
+ FOR_EACH(ichild, 0, SLN_TREE_ARITY_MAX) {
+ if(r < cumul[ichild]) {
+ leaf_proba *= proba[ichild];
+ node = sln_node_get_child(tree, node, ichild);
+ break;
+ }
+ }
+ ASSERT(ichild < SLN_TREE_ARITY_MAX); /* A node should have been sampled */
+ }
+
+exit:
+ if(out_leaf_proba) *out_leaf_proba = leaf_proba;
+ return node;
+error:
+ node = NULL;
+ leaf_proba = NaN;
+ goto exit;
+}
+
double
sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber)
{
