commit 444d1e73ceedb1ba93b60861213c8618e3bdb6e1
parent cbad96bb9d303f4897ab3acd323836459c73515c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 10 Apr 2026 10:40:09 +0200
Allow a tree arity greater than two
Currently, the maximum tree arity is set to 256.
The method for calculating the tree depth has therefore been updated, as
it was based on the assumption that the tree was binary, which is no
longer the case. Hence the update to the sln_tree_get_desc function.
The sln-build command now accepts a new (optional) argument that defines
the arity of the tree to be built. The changes made to the sln-get tool
are much more significant. Its tree traversal options were based on the
assumption that the tree was necessarily binary, which is no longer
necessarily the case. They have therefore been replaced by the -c
option, which allows the user to specify the child node to select and,
optionally, the number of levels to descend. The man pages for both
tools have been updated accordingly.
Finally, the line-sampling procedure in the sln-slab program has been
rewritten to support the tree's variable arity. At each node, a
cumulative is calculated from the list of its children. This is then
used to select a child based on its importance. Previously, this
cumulative was implicit, since the choice was between two children;
therefore, only a Bernoulli test was necessary.
Note that the performance of the calculation proposed by sln-slab can be
improved by partitioning the lines using a tree with an arity slightly
greater than 2. This behavior can be explained by the higher quality of
the polylines in each partition: they better correspond to the spectrum
they represent, since they avoid several simplification steps caused by
the successive merging of child polylines into their parent node. In
other words, a shallower tree depth implies fewer steps of polyline
merging and thus less approximation. However, despite this gain in
quality, the performance improvement is quickly offset by the increase
in arity, as the sampling procedure becomes more computationally
intensive as the number of child nodes increases (calculating a node’s
cumulative is more computationally intensive). It can therefore be
difficult to define the tree’s arity in advance to achieve the best
compromise between computational performance and memory usage.
Diffstat:
7 files changed, 191 insertions(+), 200 deletions(-)
diff --git a/doc/sln-build.1 b/doc/sln-build.1
@@ -15,7 +15,7 @@
.\"
.\" 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 April 2, 2026
+.Dd April 8, 2026
.Dt SLN-BUILD 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -26,6 +26,7 @@
.Sh SYNOPSIS
.Nm
.Op Fl hsv
+.Op Fl a Ar arity
.Op Fl e Ar polyline_opt Ns Op : Ns Ar polyline_opt No ...
.Op Fl l Ar line_profile
.Op Fl o Ar accel_struct
@@ -48,8 +49,13 @@ metadata, are formatted in HITRAN file formats.
If not provided, the list of lines is read from standard input.
.Pp
The options are as follows:
-.\""""""""""""""""""""""""""""""""""
.Bl -tag -width Ds
+.\""""""""""""""""""""""""""""""""""
+.It Fl a Ar arity
+Maximum number of children of an internal node in the constructed tree
+that partitions the input lines.
+It cannot be less than 2, which is its default value.
+.\""""""""""""""""""""""""""""""""""
.It Fl e Ar polyline_opt Ns Op : Ns Ar polyline_opt No ...
Configure the polylines, i.e., the data used to encode the shape of the
spectrum in the acceleration structure.
diff --git a/doc/sln-get.1 b/doc/sln-get.1
@@ -15,7 +15,7 @@
.\"
.\" 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 April 8, 2026
+.Dd April 10, 2026
.Dt SLN-GET 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -25,10 +25,9 @@
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh SYNOPSIS
.Nm
-.Op Fl hlmnrsv
+.Op Fl hmnsv
+.Op Fl c Ar child_id Ns Op : Ns Ar level_count
.Op Fl d Ar level
-.Op Fl L Ar nlevels
-.Op Fl R Ar nlevels
.Op Fl w Ar wavenumber
.Fl i Ar lines
.Fl p Ar molparams
@@ -58,8 +57,8 @@ or the value of the spectrum it represents
.Pq options Fl w .
By default, the node queried is the root node.
The caller can select another node by visiting the tree using the
-traversal options
-.Pq Fl L , Fl l , Fl R No and Fl r .
+traversal option
+.Pq option Fl c .
.Pp
Finally, the
.Fl d
@@ -89,22 +88,26 @@ binary, or in plain text HITRAN format, depending on whether the
.Fl s
option is set or not, respectively.
.\""""""""""""""""""""""""""""""""""
-.It Fl L Ar nlevels
-Descend the tree by visiting the left child of the nodes traversed up to
-.Ar nlevels
+.It Fl c Ar child_id Ns Op : Ns Ar level_count
+Traverse the tree downward by visiting the
+.Ar child_id Ns -th
+child of the current node.
+The value of
+.Ar child_id
+is between 0 and the number of children of the current node minus 1.
+If
+.Ar child_id
+is greater than or equal to the number of children of the node, it is
+truncated to the maximum value it could take.
+.Pp
+This traversal is repeated
+.Ar level_count
times.
-The traversal of the tree stops when the number of levels traversed is
-reached or when one of the nodes visited is a leaf.
-The node at which the traversal stopped then becomes the current node.
-.\""""""""""""""""""""""""""""""""""
-.It Fl l
-Visit the left child of the current node.
-This is a shortcut for the
-.Fl L Ns Ar 1
-option.
-The current node then becomes this left child.
-If the node has no left child, i.e., if it is a leaf, then the current
-node is not updated.
+If
+.Ar level_count
+is not provided, its default value is 1.
+The tree traversal stops when the specified number of levels has been
+reached or when one of the visited nodes is a leaf.
.\""""""""""""""""""""""""""""""""""
.It Fl m
Prints the polyline of the current node, i.e. the mesh representing the high
@@ -124,23 +127,6 @@ the number of the node's children.
Isotopologue metadata from which the tree was built.
The data is in HITRAN format.
.\""""""""""""""""""""""""""""""""""
-.It Fl R Ar nlevels
-Descend the tree by visiting the right child of the nodes traversed up
-to
-.Ar nlevels
-times.
-This option behaves in the same way as the
-.Fl L
-option, but here for the right child of the nodes, instead of their left
-child.
-.\""""""""""""""""""""""""""""""""""
-.It Fl r
-Query the right child of the current node.
-This option behaves in the same way as the
-.Fl l
-option, but here for the right child of the node, instead of the left
-child.
-.\""""""""""""""""""""""""""""""""""
.It Fl s
Specifies that input lines are formatted according to the binary format
as written by the
@@ -184,34 +170,34 @@ sln-get -i lines.par -p molparams.txt tree.sln
.Ed
.\""""""""""""""""""""""""""""""""""
.Pp
-Print the description of the child left of the tree root:
+Print the description of the first child of the tree root:
.Bd -literal -offset Ds
-sln-get -i lines.par -p molparams.txt -ln tree.sln
+sln-get -i lines.par -p molparams.txt -c0 -n tree.sln
.Ed
.\""""""""""""""""""""""""""""""""""
.Pp
-Print the description of one of the grandchildren of the root of the
-tree, in this case the right child of its left child:
+Print the description of one of the tree root's grandchildren,
+specifically the second grandchild of its first child:
.Bd -literal -offset Ds
-sln-get -i lines.par -p molparams.txt -lrn tree.sln
+sln-get -i lines.par -p molparams.txt -c0 -c1 -n tree.sln
.Ed
.\""""""""""""""""""""""""""""""""""
.Pp
-Descend the tree by first visiting the three left children of the first
+Descend the tree by first visiting the three first children of the first
three levels of the tree
-.Pq option Fl L Ns Ar 3 ,
-then the right child
-.Pq option Fl r ,
-and finally the left child
-.Pq option Fl l
-of the next two levels.
+.Pq option Fl c Ns Ar 0 : Ns Ar 3 ,
+then the second child
+.Pq option Fl c Ns Ar 1 ,
+and finally the third child
+of the next two levels
+.Pq option Fl c Ns Ar 2 : Ns Ar 2 .
Then output the polyline of the node reached
.Pq option Fl m
and save it to the
.Pa polyline.txt
file:
.Bd -literal -offset Ds
-sln-get -i lines.par -p molparams.txt -L3 -rl -m tree.sln \e
+sln-get -i lines.par -p molparams.txt -c0:3 -c1 -c2:2 -m tree.sln \e
> polyline.txt
.Ed
.\""""""""""""""""""""""""""""""""""
@@ -219,7 +205,7 @@ sln-get -i lines.par -p molparams.txt -L3 -rl -m tree.sln \e
Print the spectrum value at 50 cm^-1 for one of the great-grandchildren
of the root:
.Bd -literal -offset Ds
-sln-get -i lines.par -p molparams.txt -L3 -w 50 tree.sln
+sln-get -i lines.par -p molparams.txt -c0:3 -w 50 tree.sln
.Ed
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh SEE ALSO
diff --git a/src/sln.h b/src/sln.h
@@ -44,7 +44,7 @@
#endif
#define SLN_TREE_DEPTH_MAX 64 /* Maximum depth of a tree */
-#define SLN_TREE_ARITY_MAX 2 /* Maximum arity of a tree */
+#define SLN_TREE_ARITY_MAX 256 /* Maximum arity of a tree */
/* Forward declaration of external data structures */
struct logger;
diff --git a/src/sln_build.c b/src/sln_build.c
@@ -55,6 +55,7 @@ struct args {
/* Miscellaneous */
int quit;
int verbose;
+ unsigned arity; /* tree arity */
enum line_list_format line_format;
};
#define ARGS_DEFAULT__ { \
@@ -78,6 +79,7 @@ struct args {
/* Miscellaneous */ \
0, /* Quit */ \
0, /* Verbose */ \
+ 2, /* Tree arity */ \
LINE_LIST_HITRAN /* lines_format */ \
}
static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
@@ -103,9 +105,9 @@ static void
usage(FILE* stream)
{
fprintf(stream,
-"usage: sln-build [-hsv] [-e polyline_opt[:polyline_opt ...]] [-l line_profile]\n"
-" [-o accel_struct] -P pressure -T temperature -m molparams\n"
-" -x mixture [lines]\n");
+"usage: sln-build [-hsv] [-a arity] [-e polyline_opt[:polyline_opt ...]]\n"
+" [-l line_profile] [-o accel_struct] -P pressure -T temperature\n"
+" -m molparams -x mixture [lines]\n");
}
static res_T
@@ -223,8 +225,12 @@ args_init(struct args* args, int argc, char** argv)
*args = ARGS_DEFAULT;
- while((opt = getopt(argc, argv, "e:hl:o:P:sT:m:vx:")) != -1) {
+ while((opt = getopt(argc, argv, "a:e:hl:o:P:sT:m:vx:")) != -1) {
switch(opt) {
+ case 'a':
+ res = cstr_to_uint(optarg, &args->arity);
+ if(res == RES_OK && args->arity < 2) res = RES_BAD_ARG;
+ break;
case 'e':
res = cstr_parse_list(optarg, ':', parse_polyline_opt, args);
break;
@@ -442,6 +448,7 @@ cmd_run(struct cmd* cmd)
tree_args.nvertices_hint = cmd->args.nvertices_hint;
tree_args.mesh_decimation_err = cmd->args.mesh_decimation_err;
tree_args.mesh_type = cmd->args.mesh_type;
+ tree_args.arity = cmd->args.arity;
if(cmd->args.output) {
write_args.filename = cmd->args.output;
diff --git a/src/sln_get.c b/src/sln_get.c
@@ -46,16 +46,8 @@ struct args {
double wavenumber; /* Wave number at which the spectrum is evaluated */
- /* Step of the tree descent.
- *
- * A bit set in one of the masks determines which side to descend to move from
- * the current depth (which corresponds to the number of the set bit) to the
- * next one. Since the tree is binary, only one of the two masks is needed, as
- * the other can be deduced from the selected mask and the current depth.
- * However, calculating both descent masks allows bugs in the analysis of
- * arguments to be detected. */
- uint64_t descent_lmask; /* Left mask */
- uint64_t descent_rmask; /* Right mask */
+ /* Steps for traversing the tree */
+ unsigned descent_path[SLN_TREE_DEPTH_MAX];
unsigned depth; /* Current depth in the tree */
/* Miscellaneous */
@@ -64,7 +56,7 @@ struct args {
int verbose;
int lines_in_shtr_format;
};
-#define ARGS_DEFAULT__ {NULL,NULL,NULL,OUTPUT_TREE_DESCRIPTOR,0,0,0,0,0,0,0,0}
+#define ARGS_DEFAULT__ {NULL,NULL,NULL,OUTPUT_TREE_DESCRIPTOR,0,{0},0,0,0,0,0}
static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
struct cmd {
@@ -87,48 +79,35 @@ static void
usage(FILE* stream)
{
fprintf(stream,
-"usage: sln-get [-hlmnrsv] [-d level] [-L nlevels] [-R nlevels] [-w wavenumber]\n"
+"usage: sln-get [-hlmnrsv] [-c child_id[:level_count]] [-d level] [-w wavenumber]\n"
" -i lines -p molparams [tree]\n");
}
-static void
-tree_descent
- (struct args* args,
- const enum child child, /* side of the tree descent */
- const unsigned nlevels) /* Number of levels to descend in the tree */
+static res_T
+tree_descent(struct args* args, const char* str)
{
- uint64_t mask = 0;
- ASSERT(args);
+ unsigned path[2] = {0/*child_id*/,1/*level_count*/};
+ unsigned i=0;
+ size_t n=0;
+ res_T res = RES_OK;
+ ASSERT(args && str);
- if(nlevels == 0 || args->depth >= SLN_TREE_DEPTH_MAX)
- return; /* Nothing to descent */
+ res = cstr_to_list_uint(str, ':', path, &n, 2);
+ if(res != RES_OK) goto error;
- if(nlevels > SLN_TREE_DEPTH_MAX) {
- /* The mask can encode up to 64 levels of descent. Thus, if the number of
- * levels is greater than or equal to 64, all bits of the mask are
- * activated. */
- mask = ~0lu;
- } else {
- /* Set mask bits from bit 0 to the nlevels^th bit */
- mask = BIT_U64(nlevels)-1;
+ for(i=0; i<path[1] && args->depth<SLN_TREE_DEPTH_MAX; ++i, ++args->depth) {
+ args->descent_path[args->depth] = path[0];
}
- /* Move the mask to the current depth */
- mask <<= args->depth;
-
- /* Here we go: descent the tree */
- switch(child) {
- case LEFT: args->descent_lmask |= mask; break;
- case RIGHT: args->descent_rmask |= mask; break;
- default: FATAL("Unreachable code\n"); break;
- }
- args->depth = MMIN(args->depth + nlevels, SLN_TREE_DEPTH_MAX);
+exit:
+ return res;
+error:
+ goto exit;
}
static res_T
args_init(struct args* args, int argc, char** argv)
{
- unsigned nlvls = 0;
int opt = 0;
res_T res = RES_OK;
@@ -136,8 +115,9 @@ args_init(struct args* args, int argc, char** argv)
*args = ARGS_DEFAULT;
- while((opt = getopt(argc, argv, "d:hi:L:lmnp:R:rsvw:")) != -1) {
+ while((opt = getopt(argc, argv, "c:d:hi:mnp:svw:")) != -1) {
switch(opt) {
+ case 'c': res = tree_descent(args, optarg); break;
case 'd':
args->output_type = OUTPUT_LEVEL_DESCRIPTOR;
res = cstr_to_uint(optarg, &args->level);
@@ -147,21 +127,9 @@ args_init(struct args* args, int argc, char** argv)
args->quit = 1;
goto exit;
case 'i': args->lines = optarg; break;
- case 'L':
- if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) {
- tree_descent(args, LEFT, nlvls);
- }
- break;
- case 'l': tree_descent(args, LEFT, 1); break;
case 'm': args->output_type = OUTPUT_NODE_MESH; break;
case 'n': args->output_type = OUTPUT_NODE_DESCRIPTOR; break;
case 'p': args->molparams = optarg; break;
- case 'R':
- if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) {
- tree_descent(args, RIGHT, nlvls);
- }
- break;
- case 'r': tree_descent(args, RIGHT, 1); break;
case 's': args->lines_in_shtr_format = 1; break;
case 'v': args->verbose += (args->verbose < 3); break;
case 'w':
@@ -190,14 +158,6 @@ args_init(struct args* args, int argc, char** argv)
MANDATORY(args->lines, "line list", 'i');
#undef MANDATORY
-#ifndef NDEBUG
- {
- const uint64_t descent_mask =
- args->depth >= 64 ? ~0ul : BIT_U64(args->depth)-1;
- ASSERT((args->descent_lmask | args->descent_rmask) == descent_mask);
- }
-#endif
-
if(optind < argc) args->tree = argv[optind];
exit:
@@ -303,7 +263,7 @@ print_level_descriptor(const struct cmd* cmd)
struct {
const struct sln_node* node;
unsigned level;
- } stack[SLN_TREE_DEPTH_MAX];
+ } stack[SLN_TREE_DEPTH_MAX*(SLN_TREE_ARITY_MAX-1/*1st node's child*/)];
int istack = 0;
/* Node data */
@@ -331,12 +291,17 @@ print_level_descriptor(const struct cmd* cmd)
ASSERT(level <= cmd->args.level);
if(!sln_node_is_leaf(node) && level < cmd->args.level) {
+ const unsigned nchildren = sln_node_get_child_count(cmd->tree, node);
+ unsigned ichild = 0;
+
/* Continue down the tree */
++level;
- /* Push the right child */
- stack[istack ].node = sln_node_get_child(cmd->tree, node, 1);
- stack[istack++].level = level;
+ /* Push the node children excepted the 1st */
+ FOR_EACH(ichild, 1, nchildren) {
+ stack[istack ].node = sln_node_get_child(cmd->tree, node, ichild);
+ stack[istack++].level = level;
+ }
node = sln_node_get_child(cmd->tree, node, 0); /* Visit the left child */
@@ -374,17 +339,16 @@ get_node(const struct cmd* cmd, unsigned* node_depth/*can be NULL*/)
if(node == NULL) goto exit; /* Tree is empty */
FOR_EACH(i, 0, cmd->args.depth) {
- const uint64_t mask = BIT_U64(i);
+ unsigned nchildren = 0;
+ unsigned ichild = 0;
if(sln_node_is_leaf(node)) break;
- if(mask & cmd->args.descent_lmask) {
- node = sln_node_get_child(cmd->tree, node, 0);
- } else if(mask & cmd->args.descent_rmask) {
- node = sln_node_get_child(cmd->tree, node, 1);
- } else {
- FATAL("Unreachable code\n");
- }
+ nchildren = sln_node_get_child_count(cmd->tree, node);
+ ichild = MMIN(cmd->args.descent_path[i], nchildren-1);
+
+ node = sln_node_get_child(cmd->tree, node, ichild);
+
++depth;
}
@@ -410,7 +374,7 @@ print_node_descriptor(const struct cmd* cmd)
printf("level: %u\n", depth);
printf("#lines: %lu\n", (unsigned long)desc.nlines);
printf("#vertices: %lu\n", (unsigned long)desc.nvertices);
- printf("#nodes: %u\n", desc.nchildren);
+ printf("#children: %u\n", desc.nchildren);
exit:
return res;
@@ -518,6 +482,7 @@ print_tree_descriptor(const struct cmd* cmd)
printf("decimation error: %.4e\n", desc.mesh_decimation_err);
printf("line profile: %s\n", sln_line_profile_cstr(desc.line_profile));
printf("#lines per leaf: %lu\n", (unsigned long)desc.max_nlines_per_leaf);
+ printf("arity: %u\n", desc.arity);
exit:
return res;
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -26,6 +26,7 @@
#include <rsys/cstr.h>
#include <rsys/mem_allocator.h>
+#include <rsys/str.h>
#include <omp.h>
@@ -307,6 +308,65 @@ 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: %c%s\n",
+ ka, cumul, nu, path[0] != '\0' ? '-' : ' ', 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_line
(const struct cmd* cmd,
@@ -314,71 +374,50 @@ sample_line
const double nu/*[cm^-1]*/,
double* out_proba)
{
- char step[64] = {0}; /* For debug */
-
+ 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; ;++depth) {
- const struct sln_node* node0 = NULL;
- const struct sln_node* node1 = NULL;
- struct sln_mesh mesh = SLN_MESH_NULL;
- struct sln_mesh mesh0 = SLN_MESH_NULL;
- struct sln_mesh mesh1 = SLN_MESH_NULL;
- double proba0 = 0;
- double ka_max = 0;
- double ka_max0 = 0;
- double ka_max1 = 0;
+ for(depth=0; !sln_node_is_leaf(node); ++depth) {
double r = 0;
+ unsigned ichild = 0;
- if(sln_node_is_leaf(node)) break;
-
- node0 = sln_node_get_child(cmd->tree, node, 0);
- node1 = sln_node_get_child(cmd->tree, node, 1);
- SLN(node_get_mesh(cmd->tree, node0, &mesh0));
- SLN(node_get_mesh(cmd->tree, node1, &mesh1));
-
- ka_max0 = sln_mesh_eval(&mesh0, nu);
- ka_max1 = sln_mesh_eval(&mesh1, nu);
- /* TODO handle ka_max[0,1] == 0 */
-
- /* 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_max = sln_mesh_eval(&mesh, nu);
- if(ka_max < ka_max0 + ka_max1) {
- step[depth] = '\0';
- if(cmd->args.verbose >= 1) {
- fprintf(stderr,
- "error: ka < ka0 + ka1; %e < %e+%e; nu=%-21.20g cm^-1; node path: %c%s\n",
- ka_max, ka_max0, ka_max1, nu, step[0] != '\0' ? '-' : ' ', step);
- }
- return NULL;
- }
-
- /* Compute the probability to sample the 1st node */
- proba0 = ka_max0 / (ka_max0 + ka_max1);
+ res = build_node_cumulative
+ (cmd, node, nu, str_cget(&path), probabilities, cumulative);
+ if(res != RES_OK) goto error;
- /* Sample the child node */
+ /* Sample a child with respect to its importance */
r = ssp_rng_canonical(rng);
- if(r < proba0) {
- node = node0;
- proba *= proba0;
- step[depth] = 'l';
- } else {
- node = node1;
- proba *= (1-proba0);
- step[depth] = 'r';
+ 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;
}
static INLINE res_T
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -584,7 +584,7 @@ sln_tree_ref_put(struct sln_tree* tree)
res_T
sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc)
{
- size_t nlines_adjusted = 0;
+ const struct sln_node* node = NULL;
unsigned depth = 0;
if(!tree || !desc) return RES_BAD_ARG;
@@ -600,25 +600,13 @@ sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc)
desc->arity = tree->args.arity;
SHTR(line_list_get_size(tree->args.lines, &desc->nlines));
- nlines_adjusted = round_up_pow2(desc->nlines);
-
- for(depth=0;
- depth<SLN_TREE_DEPTH_MAX && !(BIT_U64(depth) & nlines_adjusted);
- ++depth);
-
- desc->depth = depth;
-#ifndef NDEBUG
- {
- unsigned max_depth = 0;
- const struct sln_node* node = sln_tree_get_root(tree);
- while(!sln_node_is_leaf(node)) {
- node = sln_node_get_child(tree, node, 0);
- ++max_depth;
- }
- CHK(max_depth == depth);
+ node = sln_tree_get_root(tree);
+ while(!sln_node_is_leaf(node)) {
+ node = sln_node_get_child(tree, node, 0);
+ ++depth;
}
-#endif
+ desc->depth = depth;
return RES_OK;
}
