star-line

sln_tree.c (30212B)

---

 /* 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 file is part of Star-Line.
  *
  * 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/>. */
 
 #include "sln.h"
 #include "sln_device_c.h"
 #include "sln_line.h"
 #include "sln_tree_c.h"
 
 #include <star/shtr.h>
 #include <star/ssp.h>
 
 #include <rsys/algorithm.h>
 #include <rsys/cstr.h>
 #include <rsys/math.h>
 
 #include <omp.h>
 
 struct stream {
   const char* name;
   FILE* fp;
   int intern_fp; /* Define if the stream was internally opened */
 };
 static const struct stream STREAM_NULL = {NULL, NULL, 0};
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static INLINE res_T
 check_molecule_concentration
   (const struct sln_device* sln,
    const char* caller,
    const enum shtr_molecule_id molecule_id,
    const double concentration)
 {
   ASSERT(sln && caller);
 
   if(concentration == 0) {
     /* A molecular concentration of zero is allowed, but may be a user error,
      * as 0 is the default concentration in the tree creation arguments.
      * Therefore, warn the user about this value so that they can determine
      * whether or not it is an error on their part. */
     WARN(sln, "%s: the concentration of %s is zero.\n",
       caller, shtr_molecule_cstr(molecule_id));
 
   } else if(concentration < 0) {
     /* Concentration cannot be negative... */
     ERROR(sln, "%s: invalid %s concentration: %g.\n",
       FUNC_NAME, shtr_molecule_cstr(molecule_id),
       concentration);
     return RES_BAD_ARG;
   }
 
   return RES_OK;
 }
 
 static res_T
 check_concentrations_list
   (const struct sln_device* sln,
    const char* caller,
    const double concentrations[SHTR_MAX_MOLECULE_COUNT])
 {
   double sum = 0;
   int i = 0;
   res_T res = RES_OK;
   ASSERT(sln && caller && concentrations);
 
   FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
     if(i == SHTR_MOLECULE_ID_NULL) continue;
 
     res = check_molecule_concentration(sln, caller, i, concentrations[i]);
     if(res != RES_OK) goto error;
 
     sum += concentrations[i];
   }
 
   /* The sum of molecular concentrations must be less than or equal to 1. It may
    * be less than 1 if the remaining part of the mixture is (implicitly) defined
    * as a radiatively inactive gas */
   if(sum > 1 && sum-1 > 1e-6) {
     ERROR(sln,
       "%s: the sum of molecule concentrations is greater than 1: %g\n",
       caller, sum);
     res = RES_BAD_ARG;
     goto error;
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 /* Check the consistency of the molecular concentrations */
 static INLINE res_T
 check_mixture_concentrations
   (struct sln_device* sln,
    const char* caller,
    const struct sln_tree_create_args* args)
 {
   double concentrations[SHTR_MAX_MOLECULE_COUNT] = {0};
   int i = 0;
   ASSERT(sln && caller && args);
 
   FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
     concentrations[i] = args->molecules[i].concentration;
   }
 
   return check_concentrations_list(sln, caller, concentrations);
 }
 
 /* Verify that the isotope abundance are valids */
 static res_T
 check_molecule_isotope_abundances
   (struct sln_device* sln,
    const char* caller,
    const struct sln_molecule* molecule)
 {
   int i = 0;
   double sum = 0;
   ASSERT(sln && caller && molecule);
 
   /* The isotopic abundances are the default ones. Nothing to do */
   if(!molecule->non_default_isotope_abundances) return RES_OK;
 
   /* The isotopic abundances are not the default ones.
    * Verify that they are valid ... */
   FOR_EACH(i, 0, SHTR_MAX_ISOTOPE_COUNT) {
     if(molecule->isotopes[i].abundance < 0) {
       const int isotope_id = i + 1; /* isotope id in [1, 10] */
       ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n",
         caller, isotope_id, shtr_molecule_cstr(i),
         molecule->isotopes[i].abundance);
       return RES_BAD_ARG;
     }
 
     sum += molecule->isotopes[i].abundance;
   }
 
   /* ... and that their sum equals 1 */
   if(!eq_eps(sum, 1, 1e-6)) {
     ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n",
       caller, shtr_molecule_cstr(i), sum);
     return RES_BAD_ARG;
   }
 
   return RES_OK;
 }
 
 static res_T
 check_molecules
   (struct sln_device* sln,
    const char* caller,
    const struct sln_tree_create_args* args)
 {
   char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0};
 
   size_t iline = 0;
   size_t nlines = 0;
   res_T res = RES_OK;
   ASSERT(args->lines);
 
   res = check_mixture_concentrations(sln, caller, args);
   if(res != RES_OK) goto error;
 
   /* Iterate over the lines to define which molecules has to be checked, i.e.,
    * the ones used in the mixture */
   SHTR(line_list_get_size(args->lines, &nlines));
   FOR_EACH(iline, 0, nlines) {
     struct shtr_line line = SHTR_LINE_NULL;
     const struct sln_molecule* molecule = NULL;
 
     SHTR(line_list_at(args->lines, iline, &line));
 
     /* This molecule was already checked */
     if(molecule_ok[line.molecule_id]) continue;
 
     molecule = args->molecules + line.molecule_id;
 
     if(molecule->cutoff <= 0) {
       /* ... cutoff either */
       ERROR(sln, "%s: invalid %s cutoff: %g.\n",
         caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff);
       return RES_BAD_ARG;
     }
 
     res = check_molecule_isotope_abundances(sln, caller, molecule);
     if(res != RES_OK) goto error;
 
     molecule_ok[line.molecule_id] = 1;
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static INLINE res_T
 check_pressure
   (const struct sln_device* sln,
    const char* caller,
    const double pressure /*[atm]*/)
 {
   if(pressure < 0) {
     ERROR(sln, "%s: invalid negative pressure %g atm\n", caller, pressure);
     return RES_BAD_ARG;
   }
   return RES_OK;
 }
 
 static INLINE res_T
 check_temperature
   (const struct sln_device* sln,
    const char* caller,
    const double temperature /*[K]*/)
 {
   if(temperature < 0) {
     ERROR(sln, "%s: invalid negative temperature %g K\n", caller, temperature);
     return RES_BAD_ARG;
   }
   return RES_OK;
 }
 
 static res_T
 check_sln_tree_create_args
   (struct sln_device* sln,
    const char* caller,
    const struct sln_tree_create_args* args)
 {
   res_T res = RES_OK;
   ASSERT(sln && caller);
 
   if(!args) return RES_BAD_ARG;
 
   if(!args->metadata) {
     ERROR(sln, "%s: the isotope metadata are missing.\n", caller);
     return RES_BAD_ARG;
   }
 
   if(!args->lines) {
     ERROR(sln, "%s: the list of lines is missing.\n", caller);
     return RES_BAD_ARG;
   }
 
   if((res = check_pressure(sln, caller, args->pressure)) != RES_OK) {
     return res;
   }
 
   if((res = check_temperature(sln, caller, args->temperature)) != RES_OK) {
     return res;
   }
 
   if(args->nvertices_hint == 0) {
     ERROR(sln,
       "%s: invalid hint on the number of vertices around the line center %lu.\n",
       caller, (unsigned long)args->nvertices_hint);
     return RES_BAD_ARG;
   }
 
   if(args->mesh_decimation_err < 0) {
     ERROR(sln, "%s: invalid decimation error %g.\n",
       caller, args->mesh_decimation_err);
     return RES_BAD_ARG;
   }
 
   if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) {
     ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type);
     return RES_BAD_ARG;
   }
 
   if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) {
     ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile);
     return RES_BAD_ARG;
   }
 
   if(args->arity < 2 || args->arity > SLN_TREE_ARITY_MAX) {
     ERROR(sln, "%s: invalid arity %u. It must be in [2, %d]\n",
       caller, args->arity, SLN_TREE_ARITY_MAX);
     return RES_BAD_ARG;
   }
 
   if(args->leaf_nlines < 1 || args->leaf_nlines > SLN_LEAF_NLINES_MAX) {
     ERROR(sln, "%s: invalid number of lines per leaf %u. It must be in [1, %d]\n",
       caller, args->leaf_nlines, SLN_LEAF_NLINES_MAX);
     return RES_BAD_ARG;
   }
 
   if(args->nthreads_hint == 0) {
     ERROR(sln, "%s: invalid number of threads %u\n",
       caller, args->nthreads_hint);
     return RES_BAD_ARG;
   }
 
   res = check_molecules(sln, caller, args);
   if(res != RES_OK) return res;
 
   return RES_OK;
 }
 
 static res_T
 check_sln_tree_read_args
   (struct sln_device* sln,
    const char* caller,
    const struct sln_tree_read_args* args)
 {
   if(!args) return RES_BAD_ARG;
 
   if(!args->metadata) {
     ERROR(sln, "%s: the isotope metadata are missing.\n", caller);
     return RES_BAD_ARG;
   }
 
   if(!args->lines) {
     ERROR(sln, "%s: the list of lines is missing.\n", caller);
     return RES_BAD_ARG;
   }
 
   if(!args->file && !args->filename) {
     ERROR(sln,
       "%s: the source file is missing. No file name or stream is provided.\n",
       caller);
     return RES_BAD_ARG;
   }
 
   return RES_OK;
 }
 
 static res_T
 check_sln_tree_write_args
   (struct sln_device* sln,
    const char* caller,
    const struct sln_tree_write_args* args)
 {
   if(!args) return RES_BAD_ARG;
 
   if(!args->file && !args->filename) {
     ERROR(sln,
       "%s: the destination file is missing. "
       "No file name or stream is provided.\n",
       caller);
     return RES_BAD_ARG;
   }
 
   return RES_OK;
 }
 
 static res_T
 check_line_thermo_props
   (const struct sln_tree* tree,
    const char* caller,
    const size_t iline,
    const struct sln_thermo_props* props)
 {
   struct shtr_line line = SHTR_LINE_NULL;
   res_T res = RES_OK;
   ASSERT(tree && caller);
 
   if(!props) goto exit; /* Default thermo props */
 
   SHTR(line_list_at(tree->args.lines, iline, &line));
 
   res = check_molecule_concentration(tree->sln, caller, line.molecule_id,
     props->concentrations[line.molecule_id]);
   if(res != RES_OK) goto error;
 
   res = check_pressure(tree->sln, caller, props->pressure);
   if(res != RES_OK) goto error;
 
   res = check_temperature(tree->sln, caller, props->temperature);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static INLINE res_T
 check_sln_thermo_props
   (const struct sln_device* sln,
    const char* caller,
    const struct sln_thermo_props* props)
 {
   res_T res = RES_OK;
   ASSERT(sln && caller);
 
   if(!props) goto exit; /* Default thermo props */
 
   res = check_concentrations_list(sln, caller, props->concentrations);
   if(res != RES_OK) goto error;
 
   res = check_pressure(sln, caller, props->pressure);
   if(res != RES_OK) goto error;
 
   res = check_temperature(sln, caller, props->temperature);
   if(res != RES_OK) goto error;
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 static INLINE void
 stream_release(struct stream* stream)
 {
   ASSERT(stream);
   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
 }
 
 static res_T
 stream_init
   (struct sln_device* sln,
    const char* caller,
    const char* name, /* NULL <=> default stream name */
    FILE* fp, /* NULL <=> open file "name" */
    const char* mode, /* mode in fopen */
    struct stream* stream)
 {
   res_T res = RES_OK;
 
   ASSERT(sln && caller && stream);
   ASSERT(fp || (name && mode));
 
   *stream = STREAM_NULL;
 
   if(fp) {
     stream->intern_fp = 0;
     stream->name = name ? name : "stream";
     stream->fp = fp;
 
   } else {
     stream->intern_fp = 1;
     stream->name = name;
     if(!(stream->fp = fopen(name, mode))) {
       ERROR(sln, "%s:%s: error opening file -- %s\n",
         caller, name, strerror(errno));
       res = RES_IO_ERR;
       goto error;
     }
   }
 
 exit:
   return res;
 error:
   if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0);
   goto exit;
 }
 
 static res_T
 create_tree
   (struct sln_device* sln,
    const char* caller,
    struct sln_tree** out_tree)
 {
   struct sln_tree* tree = NULL;
   res_T res = RES_OK;
   ASSERT(sln && caller && out_tree);
 
   tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree));
   if(!tree) {
     ERROR(sln, "%s: could not allocate the tree data structure.\n",
       caller);
     res = RES_MEM_ERR;
     goto error;
   }
   ref_init(&tree->ref);
   SLN(device_ref_get(sln));
   tree->sln = sln;
   darray_node_init(sln->allocator, &tree->nodes);
   darray_vertex_init(sln->allocator, &tree->vertices);
 
 exit:
   *out_tree = tree;
   return res;
 error:
   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
   goto exit;
 }
 
 static INLINE int
 cmp_nu_vtx(const void* key, const void* item)
 {
   const float nu = *((const float*)key);
   const struct sln_vertex* vtx = item;
   if(nu < vtx->wavenumber) return -1;
   if(nu > vtx->wavenumber) return +1;
   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)
 {
   struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref);
   struct sln_device* sln = NULL;
   ASSERT(ref);
   sln = tree->sln;
   darray_node_release(&tree->nodes);
   darray_vertex_release(&tree->vertices);
   if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines));
   if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata));
   MEM_RM(sln->allocator, tree);
   SLN(device_ref_put(sln));
 }
 
 /*******************************************************************************
  * Local function
  ******************************************************************************/
 unsigned
 node_child_count(const struct sln_node* node, const unsigned tree_arity)
 {
   size_t nlines = 0; /* #lines in the node */
   size_t nlines_per_child =  0; /* Max #lines in a child */
   size_t nchildren = 0;
 
   /* Pre-conditions */
   ASSERT(node && tree_arity >= 2);
 
   /* Retrieve the node data and compute the #lines it partitions */
   nlines = node->range[1] - node->range[0] + 1;
   ASSERT(nlines);
 
   /* Based on the arity of the tree, calculate how the lines of the node are
    * distributed among its children. For low lines count, i.e. when the minimum
    * number of lines par child is less than the tree arity, the policy below
    * prioritizes an equal distribution of lines among the children over
    * maintaining the tree's arity.  Thus, if a smaller number of children
    * results in a more equitable distribution, this option is preferred over
    * ensuring a number of children equal to the tree's arity. In other words,
    * the tree's balance is prioritized. */
   nlines_per_child = (nlines + tree_arity-1/*ceil*/)/tree_arity;
 
   /* From the previous line repartition, compute the number of children */
   nchildren = (nlines + nlines_per_child-1/*ceil*/)/nlines_per_child;
   ASSERT(nchildren >= 2);
 
   ASSERT(nchildren <= UINT_MAX);
   return (unsigned)nchildren;
 }
 
 /*******************************************************************************
  * Exported symbols
  ******************************************************************************/
 res_T
 sln_tree_create
   (struct sln_device* device,
    const struct sln_tree_create_args* args,
    struct sln_tree** out_tree)
 {
   struct sln_tree* tree = NULL;
   unsigned nthreads_max = 0;
   res_T res = RES_OK;
 
   if(!device || !out_tree) { res = RES_BAD_ARG; goto error; }
   res = check_sln_tree_create_args(device, FUNC_NAME, args);
   if(res != RES_OK) goto error;
 
   res = create_tree(device, FUNC_NAME, &tree);
   if(res != RES_OK) goto error;
   SHTR(line_list_ref_get(args->lines));
   SHTR(isotope_metadata_ref_get(args->metadata));
   tree->args = *args;
 
   /* Set the #threads to match the maximum number of available threads */
   nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs());
   tree->args.nthreads_hint = MMIN(tree->args.nthreads_hint, nthreads_max);
 
   res = tree_build(tree);
   if(res != RES_OK) goto error;
 
 exit:
   if(out_tree) *out_tree = tree;
   return res;
 error:
   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
   goto exit;
 }
 
 res_T
 sln_tree_read
   (struct sln_device* sln,
    const struct sln_tree_read_args* args,
    struct sln_tree** out_tree)
 {
   hash256_T hash_mdata1;
   hash256_T hash_mdata2;
   hash256_T hash_lines1;
   hash256_T hash_lines2;
 
   struct stream stream = STREAM_NULL;
   struct sln_tree* tree = NULL;
   size_t n = 0;
   int version = 0;
   res_T res = RES_OK;
 
   if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; }
   res = check_sln_tree_read_args(sln, FUNC_NAME, args);
   if(res != RES_OK) goto error;
 
   res = create_tree(sln, FUNC_NAME, &tree);
   if(res != RES_OK) goto error;
 
   res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream);
   if(res != RES_OK) goto error;
 
   #define READ(Var, Nb) {                                                      \
     if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {                \
       if(feof(stream.fp)) {                                                    \
         res = RES_BAD_ARG;                                                     \
       } else if(ferror(stream.fp)) {                                           \
         res = RES_IO_ERR;                                                      \
       } else {                                                                 \
         res = RES_UNKNOWN_ERR;                                                 \
       }                                                                        \
       ERROR(sln, "%s: error loading the tree structure -- %s.\n",              \
         stream.name, res_to_cstr(res));                                        \
       goto error;                                                              \
     }                                                                          \
   } (void)0
   READ(&version, 1);
   if(version != SLN_TREE_VERSION) {
     ERROR(sln,
       "%s: unexpected tree version %d. Expecting a tree in version %d.\n",
       stream.name, version, SLN_TREE_VERSION);
     res = RES_BAD_ARG;
     goto error;
   }
 
   res = shtr_isotope_metadata_hash(args->metadata, hash_mdata1);
   if(res != RES_OK) goto error;
 
   READ(hash_mdata2, sizeof(hash256_T));
   if(!hash256_eq(hash_mdata1, hash_mdata2)) {
     ERROR(sln,
       "%s: the input isotopic metadata are not those used "
       "during tree construction.\n", stream.name);
     res = RES_BAD_ARG;
     goto error;
   }
 
   SHTR(isotope_metadata_ref_get(args->metadata));
   tree->args.metadata = args->metadata;
 
   READ(hash_lines1, sizeof(hash256_T));
   if(!args->disable_line_hash_check) {
     res = shtr_line_list_hash(args->lines, hash_lines2);
     if(res != RES_OK) goto error;
 
     if(!hash256_eq(hash_lines1, hash_lines2)) {
       ERROR(sln,
         "%s: the input list of lines is not the one used to build the tree.\n",
         stream.name);
       res = RES_BAD_ARG;
       goto error;
     }
   }
 
   SHTR(line_list_ref_get(args->lines));
   tree->args.lines = args->lines;
 
   READ(&n, 1);
   if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error;
   READ(darray_node_data_get(&tree->nodes), n);
 
   READ(&n, 1);
   if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error;
   READ(darray_vertex_data_get(&tree->vertices), n);
 
   READ(&tree->args.line_profile, 1);
   READ(&tree->args.molecules, 1);
   READ(&tree->args.pressure, 1);
   READ(&tree->args.temperature, 1);
   READ(&tree->args.nvertices_hint, 1);
   READ(&tree->args.mesh_decimation_err, 1);
   READ(&tree->args.mesh_type, 1);
   READ(&tree->args.arity, 1);
   READ(&tree->args.leaf_nlines, 1);
   #undef READ
 
 exit:
   stream_release(&stream);
   if(out_tree) *out_tree = tree;
   return res;
 error:
   if(tree) { SLN(tree_ref_put(tree)); tree = NULL; }
   goto exit;
 }
 
 res_T
 sln_tree_ref_get(struct sln_tree* tree)
 {
   if(!tree) return RES_BAD_ARG;
   ref_get(&tree->ref);
   return RES_OK;
 }
 
 res_T
 sln_tree_ref_put(struct sln_tree* tree)
 {
   if(!tree) return RES_BAD_ARG;
   ref_put(&tree->ref, release_tree);
   return RES_OK;
 }
 
 res_T
 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc)
 {
   const struct sln_node* node = NULL;
   unsigned depth = 0;
 
   if(!tree || !desc) return RES_BAD_ARG;
 
   desc->mesh_decimation_err = tree->args.mesh_decimation_err;
   desc->mesh_type = tree->args.mesh_type;
   desc->line_profile = tree->args.line_profile;
   desc->nnodes = darray_node_size_get(&tree->nodes);
   desc->nvertices = darray_vertex_size_get(&tree->vertices);
   desc->pressure = tree->args.pressure; /* [atm] */
   desc->temperature = tree->args.temperature; /* [K] */
   desc->arity = tree->args.arity;
   desc->leaf_nlines = tree->args.leaf_nlines;
 
   SHTR(line_list_get_size(tree->args.lines, &desc->nlines));
 
   node = sln_tree_get_root(tree);
   while(!sln_node_is_leaf(node)) {
     node = sln_node_get_child(tree, node, 0);
     ++depth;
   }
   desc->depth = depth;
 
   return RES_OK;
 }
 
 const struct sln_node*
 sln_tree_get_root(const struct sln_tree* tree)
 {
   ASSERT(tree);
   if(darray_node_size_get(&tree->nodes)) {
     return darray_node_cdata_get(&tree->nodes);
   } else {
     return NULL;
   }
 }
 
 res_T
 sln_tree_get_line
   (const struct sln_tree* tree,
    const size_t iline,
    const struct sln_thermo_props* props,
    struct sln_line* line)
 {
   size_t nlines = 0;
   res_T res = RES_OK;
 
   if(!tree || !line) { res = RES_BAD_ARG; goto error; }
 
   SHTR(line_list_get_size(tree->args.lines, &nlines));
   if(iline >= nlines) { res = RES_BAD_ARG; goto error; }
 
   res = check_line_thermo_props(tree, FUNC_NAME, iline, props);
   if(res != RES_OK) goto error;
 
   res = line_setup(tree, iline, props, line);
   if(res != RES_OK) {
     ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n",
       FUNC_NAME, iline, res_to_cstr(res));
     goto error;
   }
 
 exit:
   return res;
 error:
   goto exit;
 }
 
 int
 sln_node_is_leaf(const struct sln_node* node)
 {
   ASSERT(node);
   return node->offset == 0;
 }
 
 unsigned
 sln_node_get_child_count
   (const struct sln_tree* tree,
    const struct sln_node* node)
 {
   ASSERT(tree && node);
 
   if(sln_node_is_leaf(node)) {
     return 0; /* No child */
   } else {
     return node_child_count(node, tree->args.arity);
   }
 }
 
 const struct sln_node*
 sln_node_get_child
   (const struct sln_tree* tree,
    const struct sln_node* node,
    const unsigned ichild)
 {
   ASSERT(node && ichild < sln_node_get_child_count(tree, node));
   ASSERT(!sln_node_is_leaf(node));
   (void)tree;
   return node + node->offset + ichild;
 }
 
 res_T
 sln_node_get_mesh
   (const struct sln_tree* tree,
    const struct sln_node* node,
    struct sln_mesh* mesh)
 {
   if(!tree || !node || !mesh) return RES_BAD_ARG;
   mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex;
   mesh->nvertices = node->nvertices;
   return RES_OK;
 }
 
 double
 sln_node_eval
   (const struct sln_tree* tree,
    const struct sln_node* node,
    const struct sln_thermo_props* props,
    const double nu)
 {
   double ka = 0;
   size_t iline;
   ASSERT(tree && node);
   ASSERT(check_sln_thermo_props(tree->sln, FUNC_NAME, props) == RES_OK);
 
   FOR_EACH(iline, node->range[0], node->range[1]+1) {
     struct sln_line line = SLN_LINE_NULL;
     res_T res = RES_OK;
 
     res = line_setup(tree, iline, props, &line);
     if(res != RES_OK) {
       WARN(tree->sln, "%s: could not setup the line %lu-- %s\n",
         FUNC_NAME, iline, res_to_cstr(res));
       continue;
     }
 
     ka += sln_line_eval(tree, &line, nu);
   }
   return ka;
 }
 
 res_T
 sln_node_get_desc
   (const struct sln_tree* tree,
    const struct sln_node* node,
    struct sln_node_desc* desc)
 {
   if(!tree || !node || !desc) return RES_BAD_ARG;
   desc->ilines[0] = node->range[0];
   desc->ilines[1] = node->range[1];
   desc->nvertices = node->nvertices;
   desc->nchildren = sln_node_get_child_count(tree, node);
   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)
 {
   const struct sln_vertex* vtx0 = NULL;
   const struct sln_vertex* vtx1 = NULL;
   const float nu = (float)wavenumber;
   size_t n; /* #vertices */
   double u; /* Linear interpolation parameter */
   ASSERT(mesh && mesh->nvertices);
 
   n = mesh->nvertices;
 
   /* Handle special cases */
   if(n == 1) return mesh->vertices[0].ka;
   if(nu < mesh->vertices[0].wavenumber
   || nu > mesh->vertices[n-1].wavenumber) {
     return 0;
   }
   if(nu == mesh->vertices[0].wavenumber) return mesh->vertices[0].ka;
   if(nu == mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka;
 
   /* Dichotomic search of the mesh vertex whose wavenumber is greater than or
    * equal to the submitted wavenumber 'nu' */
   vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx);
   vtx0 = vtx1 - 1;
   ASSERT(vtx1); /* A vertex is necessary found ...*/
   ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */
   ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber);
 
   /* Compute the linear interpolation parameter */
   u = (wavenumber - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber);
   u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */
 
   if(u == 0) return vtx0->ka;
   if(u == 1) return vtx1->ka;
   return u*(vtx1->ka - vtx0->ka) + vtx0->ka;
 }
 
 res_T
 sln_tree_write
   (const struct sln_tree* tree,
    const struct sln_tree_write_args* args)
 {
   struct stream stream = STREAM_NULL;
   size_t nnodes, nverts;
   hash256_T hash_mdata;
   hash256_T hash_lines;
   res_T res = RES_OK;
 
   if(!tree) { res = RES_BAD_ARG; goto error; }
   res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args);
   if(res != RES_OK) goto error;
 
   res = shtr_isotope_metadata_hash(tree->args.metadata, hash_mdata);
   if(res != RES_OK) goto error;
   res = shtr_line_list_hash(tree->args.lines, hash_lines);
   if(res != RES_OK) goto error;
 
   res = stream_init
     (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream);
   if(res != RES_OK) goto error;
 
   #define WRITE(Var, Nb) {                                                     \
     if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) {               \
       ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n",                \
         FUNC_NAME, stream.name, strerror(errno));                              \
       res = RES_IO_ERR;                                                        \
       goto error;                                                              \
     }                                                                          \
   } (void)0
   WRITE(&SLN_TREE_VERSION, 1);
   WRITE(hash_mdata, sizeof(hash256_T));
   WRITE(hash_lines, sizeof(hash256_T));
 
   nnodes = darray_node_size_get(&tree->nodes);
   WRITE(&nnodes, 1);
   WRITE(darray_node_cdata_get(&tree->nodes), nnodes);
 
   nverts = darray_vertex_size_get(&tree->vertices);
   WRITE(&nverts, 1);
   WRITE(darray_vertex_cdata_get(&tree->vertices), nverts);
 
   WRITE(&tree->args.line_profile, 1);
   WRITE(&tree->args.molecules, 1);
   WRITE(&tree->args.pressure, 1);
   WRITE(&tree->args.temperature, 1);
   WRITE(&tree->args.nvertices_hint, 1);
   WRITE(&tree->args.mesh_decimation_err, 1);
   WRITE(&tree->args.mesh_type, 1);
   WRITE(&tree->args.arity, 1);
   WRITE(&tree->args.leaf_nlines, 1);
   #undef WRITE
 
 exit:
   stream_release(&stream);
   return res;
 error:
   goto exit;
 }