star-line

Structure for accelerating line importance sampling
git clone git://git.meso-star.fr/star-line.git
Log | Files | Refs | README | LICENSE

test_sln_tree.c (16191B)


      1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2026 Université de Lorraine
      3  * Copyright (C) 2022 Centre National de la Recherche Scientifique
      4  * Copyright (C) 2022 Université Paul Sabatier
      5  *
      6  * This program is free software: you can redistribute it and/or modify
      7  * it under the terms of the GNU General Public License as published by
      8  * the Free Software Foundation, either version 3 of the License, or
      9  * (at your option) any later version.
     10  *
     11  * This program is distributed in the hope that it will be useful,
     12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     14  * GNU General Public License for more details.
     15  *
     16  * You should have received a copy of the GNU General Public License
     17  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     18 
     19 #include "test_sln_lines.h"
     20 #include "sln.h"
     21 
     22 #include <star/shtr.h>
     23 
     24 #include <rsys/algorithm.h>
     25 #include <rsys/math.h>
     26 #include <rsys/mem_allocator.h>
     27 
     28 /*******************************************************************************
     29  * Helper function
     30  ******************************************************************************/
     31 /* Return the index of the line in the line_list or SIZE_MAX if the line does
     32  * not exist */
     33 static INLINE size_t
     34 find_line
     35   (struct shtr_line_list* line_list,
     36    const struct shtr_line* line)
     37 {
     38   struct shtr_line ln = SHTR_LINE_NULL;
     39   size_t lo, hi, mid;
     40   size_t iline;
     41   size_t nlines;
     42 
     43   CHK(shtr_line_list_get_size(line_list, &nlines) == RES_OK);
     44 
     45   /* Dichotomic search */
     46   lo = 0; hi = nlines -1;
     47   while(lo < hi) {
     48     mid = (lo+hi)/2;
     49 
     50     CHK(shtr_line_list_at(line_list, mid, &ln) == RES_OK);
     51     if(line->wavenumber > ln.wavenumber) {
     52       lo = mid + 1;
     53     } else {
     54       hi = mid;
     55     }
     56   }
     57   iline = lo;
     58 
     59   CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK);
     60   if(ln.wavenumber != line->wavenumber) return SIZE_MAX;
     61 
     62 
     63   /* Find a line with the same wavenumber as the one searched for and whose
     64    * other member variables are also equal to those of the line searched for */
     65   while(ln.wavenumber == line->wavenumber
     66      && !shtr_line_eq(&ln, line)
     67      && iline < nlines) {
     68     iline += 1;
     69     CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK);
     70   }
     71 
     72   return shtr_line_eq(&ln, line) ? iline : SIZE_MAX;
     73 }
     74 
     75 /* This test assumes that all the lines contained into the list are
     76  * partitionned in the tree */
     77 static void
     78 check_tree_lines
     79   (const struct sln_tree* tree,
     80    struct shtr_line_list* line_list)
     81 {
     82   #define STACK_SIZE \
     83     ( SLN_TREE_DEPTH_MAX \
     84     * SLN_TREE_ARITY_MAX - 1/*1st child*/ \
     85     + 1/*Dummy node*/)
     86   const struct sln_node* stack[STACK_SIZE] = {NULL};
     87   struct sln_tree_desc desc = SLN_TREE_DESC_NULL;
     88   size_t istack = 0;
     89   const struct sln_node* node = NULL;
     90 
     91   char* found_lines = NULL;
     92   size_t found_nlines = 0;
     93   size_t line_list_sz;
     94 
     95   CHK(shtr_line_list_get_size(line_list, &line_list_sz) == RES_OK);
     96   CHK(sln_tree_get_desc(tree, &desc) == RES_OK);
     97 
     98   CHK(found_lines = mem_calloc(line_list_sz, sizeof(char)));
     99 
    100   stack[istack++] = NULL; /* Dummy node to stop the recursion */
    101 
    102   node = sln_tree_get_root(tree);
    103   while(node) {
    104     if(!sln_node_is_leaf(node)) {
    105       unsigned i = 0;
    106       unsigned n = sln_node_get_child_count(tree, node);
    107 
    108       CHK(sln_node_get_line_count(node) > desc.max_nlines_per_leaf);
    109 
    110       FOR_EACH(i, 1, n) {
    111         stack[istack++] = sln_node_get_child(tree, node, i);
    112       }
    113       node = sln_node_get_child(tree, node, 0); /* Handle the child 0 */
    114 
    115     } else {
    116       size_t iline, nlines;
    117 
    118       nlines = sln_node_get_line_count(node);
    119       CHK(nlines <= desc.max_nlines_per_leaf);
    120       FOR_EACH(iline, 0, nlines) {
    121         struct shtr_line line = SHTR_LINE_NULL;
    122         size_t found_iline = 0;
    123         CHK(sln_node_get_line(tree, node, iline, &line) == RES_OK);
    124         found_iline = find_line(line_list, &line);
    125         CHK(found_iline != SIZE_MAX); /* Line is found */
    126 
    127         if(!found_lines[found_iline]) {
    128           found_nlines += 1;
    129           found_lines[found_iline] = 1;
    130         }
    131       }
    132 
    133       node = stack[--istack]; /* Pop the next node */
    134     }
    135   }
    136 
    137   /* Check that all lines are found */
    138   CHK(found_nlines == line_list_sz);
    139 
    140   mem_rm(found_lines);
    141   #undef STACK_SIZE
    142 }
    143 
    144 static void
    145 test_tree
    146   (struct sln_device* sln,
    147    const struct sln_tree_create_args* tree_args_in,
    148    struct shtr_line_list* line_list)
    149 {
    150   struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT;
    151   struct sln_tree_desc desc = SLN_TREE_DESC_NULL;
    152   struct sln_node_desc node_desc = SLN_NODE_DESC_NULL;
    153   struct sln_mesh mesh = SLN_MESH_NULL;
    154   struct sln_tree* tree = NULL;
    155   const struct sln_node* node = NULL;
    156   struct shtr_line line = SHTR_LINE_NULL;
    157   size_t nlines = 0;
    158 
    159   CHK(sln && tree_args_in && line_list);
    160   tree_args = *tree_args_in;
    161 
    162   CHK(sln_tree_create(NULL, &tree_args, &tree) == RES_BAD_ARG);
    163   CHK(sln_tree_create(sln, NULL, &tree) == RES_BAD_ARG);
    164   CHK(sln_tree_create(sln, &tree_args, NULL) == RES_BAD_ARG);
    165   CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK);
    166 
    167   CHK(shtr_line_list_get_size(line_list, &nlines) == RES_OK);
    168 
    169   CHK(sln_tree_get_desc(NULL, &desc) == RES_BAD_ARG);
    170   CHK(sln_tree_get_desc(tree, NULL) == RES_BAD_ARG);
    171   CHK(sln_tree_get_desc(tree, &desc) == RES_OK);
    172 
    173   CHK(desc.max_nlines_per_leaf >= 1);
    174   CHK(desc.mesh_decimation_err == tree_args.mesh_decimation_err);
    175   CHK(desc.mesh_type == tree_args.mesh_type);
    176   CHK(desc.line_profile == tree_args.line_profile);
    177   CHK(desc.nlines == nlines);
    178   CHK(desc.nvertices >= 1);
    179   CHK(desc.nnodes >= 1);
    180   CHK(desc.depth
    181   == (unsigned)ceil(log((double)desc.nlines)/log((double)desc.arity)));
    182   CHK(desc.pressure == tree_args.pressure);
    183   CHK(desc.temperature == tree_args.temperature);
    184   CHK(desc.arity == tree_args.arity);
    185 
    186   CHK(node = sln_tree_get_root(tree));
    187   CHK(node != NULL);
    188 
    189   CHK(sln_node_get_desc(NULL, node, &node_desc) == RES_BAD_ARG);
    190   CHK(sln_node_get_desc(tree, NULL, &node_desc) == RES_BAD_ARG);
    191   CHK(sln_node_get_desc(tree, node, NULL) == RES_BAD_ARG);
    192   CHK(sln_node_get_desc(tree, node, &node_desc) == RES_OK);
    193   CHK(node_desc.nlines = desc.nlines);
    194   CHK(node_desc.nvertices >= 1 && node_desc.nvertices < desc.nvertices);
    195   CHK(node_desc.nchildren <= desc.arity);
    196 
    197   while(!sln_node_is_leaf(node)) {
    198     struct sln_node_desc node_desc_next = SLN_NODE_DESC_NULL;
    199 
    200     CHK(sln_node_get_child_count(tree, node) != 0);
    201     CHK(sln_node_get_child_count(tree, node) <= desc.arity);
    202 
    203     node = sln_node_get_child(tree, node, 0);
    204 
    205     CHK(sln_node_get_desc(tree, node, &node_desc_next) == RES_OK);
    206     CHK(node_desc_next.nlines >= 1);
    207     CHK(node_desc_next.nlines < node_desc.nlines);
    208     CHK(node_desc_next.nvertices >= 1);
    209     if(sln_node_is_leaf(node)) {
    210       CHK(node_desc_next.nchildren == 0);
    211     } else {
    212       CHK(node_desc_next.nchildren != 0);
    213       CHK(node_desc_next.nchildren <= desc.arity);
    214     }
    215 
    216     node_desc = node_desc_next;
    217   }
    218   CHK(sln_node_get_child_count(tree, node) == 0);
    219 
    220   CHK(sln_node_get_line_count(node) <= desc.max_nlines_per_leaf);
    221   CHK(sln_node_get_line(NULL, node, 0, &line) == RES_BAD_ARG);
    222   CHK(sln_node_get_line(tree, NULL, 0, &line) == RES_BAD_ARG);
    223   CHK(sln_node_get_line(tree, node, desc.max_nlines_per_leaf+1, &line)
    224     == RES_BAD_ARG);
    225   CHK(sln_node_get_line(tree, node, 0, NULL) == RES_BAD_ARG);
    226   CHK(sln_node_get_line(tree, node, 0, &line) == RES_OK);
    227   CHK(sln_node_get_line(tree, sln_tree_get_root(tree), 0, &line) == RES_OK);
    228 
    229   CHK(sln_node_get_mesh(NULL, node, &mesh) == RES_BAD_ARG);
    230   CHK(sln_node_get_mesh(tree, NULL, &mesh) == RES_BAD_ARG);
    231   CHK(sln_node_get_mesh(tree, node, NULL) == RES_BAD_ARG);
    232   CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK);
    233 
    234   CHK(node = sln_tree_get_root(tree));
    235   check_tree_lines(tree, line_list);
    236 
    237   CHK(sln_tree_ref_get(NULL) == RES_BAD_ARG);
    238   CHK(sln_tree_ref_get(tree) == RES_OK);
    239   CHK(sln_tree_ref_put(NULL) == RES_BAD_ARG);
    240   CHK(sln_tree_ref_put(tree) == RES_OK);
    241   CHK(sln_tree_ref_put(tree) == RES_OK);
    242 
    243   tree_args.mesh_decimation_err = -1;
    244   CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG);
    245 
    246   tree_args.mesh_decimation_err = 1e-1f;
    247   tree_args.mesh_type = SLN_MESH_TYPES_COUNT__;
    248   CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG);
    249 
    250   tree_args.mesh_type = SLN_MESH_FIT;
    251   tree_args.line_profile = SLN_LINE_PROFILES_COUNT__;
    252   CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG);
    253 
    254   tree_args.line_profile = SLN_LINE_PROFILE_VOIGT;
    255   CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK);
    256   CHK(sln_tree_ref_put(tree) == RES_OK);
    257 }
    258 
    259 static void
    260 check_mesh_equality(const struct sln_mesh* mesh1, const struct sln_mesh* mesh2)
    261 {
    262   size_t i;
    263   CHK(mesh1 && mesh2);
    264   CHK(mesh1->nvertices == mesh2->nvertices);
    265 
    266   FOR_EACH(i, 0, mesh1->nvertices) {
    267     CHK(mesh1->vertices[i].wavenumber == mesh2->vertices[i].wavenumber);
    268     CHK(mesh1->vertices[i].ka == mesh2->vertices[i].ka);
    269   }
    270 }
    271 
    272 static void
    273 check_node_equality
    274   (const struct sln_tree* tree1,
    275    const struct sln_tree* tree2,
    276    const struct sln_node* node1,
    277    const struct sln_node* node2)
    278 {
    279   struct sln_mesh mesh1 = SLN_MESH_NULL;
    280   struct sln_mesh mesh2 = SLN_MESH_NULL;
    281   size_t iline, nlines;
    282 
    283   CHK(node1 && node2);
    284   CHK(sln_node_is_leaf(node1) == sln_node_is_leaf(node2));
    285   CHK(sln_node_get_line_count(node1) == sln_node_get_line_count(node2));
    286   nlines = sln_node_get_line_count(node1);
    287 
    288   FOR_EACH(iline, 0, nlines) {
    289     struct shtr_line line1 = SHTR_LINE_NULL;
    290     struct shtr_line line2 = SHTR_LINE_NULL;
    291     CHK(sln_node_get_line(tree1, node1, iline, &line1) == RES_OK);
    292     CHK(sln_node_get_line(tree2, node2, iline, &line2) == RES_OK);
    293     CHK(shtr_line_eq(&line1, &line2));
    294   }
    295 
    296   CHK(sln_node_get_mesh(tree1, node1, &mesh1) == RES_OK);
    297   CHK(sln_node_get_mesh(tree2, node2, &mesh2) == RES_OK);
    298   check_mesh_equality(&mesh1, &mesh2);
    299 }
    300 
    301 static void
    302 check_tree_equality
    303   (const struct sln_tree* tree1,
    304    const struct sln_tree* tree2)
    305 {
    306   #define STACK_SIZE \
    307     (  SLN_TREE_DEPTH_MAX \
    308     * (SLN_TREE_ARITY_MAX-1/*1st child*/) * 2/*#trees*/ \
    309     + 1/*Dummy node*/)
    310 
    311   const struct sln_node* stack[STACK_SIZE] = {NULL};
    312   int istack = 0;
    313 
    314   struct sln_tree_desc desc1 = SLN_TREE_DESC_NULL;
    315   struct sln_tree_desc desc2 = SLN_TREE_DESC_NULL;
    316   const struct sln_node* node1 = NULL;
    317   const struct sln_node* node2 = NULL;
    318 
    319   CHK(sln_tree_get_desc(tree1, &desc1) == RES_OK);
    320   CHK(sln_tree_get_desc(tree2, &desc2) == RES_OK);
    321   CHK(desc1.max_nlines_per_leaf == desc2.max_nlines_per_leaf);
    322   CHK(desc1.mesh_decimation_err == desc2.mesh_decimation_err);
    323   CHK(desc1.mesh_type == desc2.mesh_type);
    324   CHK(desc1.line_profile == desc2.line_profile);
    325   CHK(desc1.pressure == desc2.pressure);
    326   CHK(desc1.temperature == desc2.temperature);
    327   CHK(desc1.depth == desc2.depth);
    328   CHK(desc1.nlines == desc2.nlines);
    329   CHK(desc1.nvertices == desc2.nvertices);
    330   CHK(desc1.nnodes == desc2.nnodes);
    331   CHK(desc1.arity == desc2.arity);
    332 
    333   stack[istack++] = NULL;
    334   stack[istack++] = NULL;
    335 
    336   node1 = sln_tree_get_root(tree1);
    337   node2 = sln_tree_get_root(tree2);
    338 
    339   while(node1 || node2) {
    340     CHK((!node1 && !node2) || (node1 && node2));
    341     check_node_equality(tree1, tree2, node1, node2);
    342 
    343     if(sln_node_is_leaf(node1)) {
    344       node2 = stack[--istack];
    345       node1 = stack[--istack];
    346     } else {
    347       unsigned i = 0;
    348       unsigned n = sln_node_get_child_count(tree1, node1);
    349 
    350       FOR_EACH(i, 1, n) {
    351         stack[istack++] = sln_node_get_child(tree1, node1, i);
    352         stack[istack++] = sln_node_get_child(tree2, node2, i);
    353       }
    354       node1 = sln_node_get_child(tree1, node1, 0);
    355       node2 = sln_node_get_child(tree2, node2, 0);
    356     }
    357   }
    358 
    359   #undef STACK_SIZE
    360 }
    361 
    362 static void
    363 test_tree_serialization
    364   (struct sln_device* sln,
    365    const struct sln_tree_create_args* tree_args)
    366 {
    367   struct sln_tree_write_args wargs = SLN_TREE_WRITE_ARGS_NULL;
    368   struct sln_tree_read_args rargs = SLN_TREE_READ_ARGS_NULL;
    369   struct sln_tree* tree1 = NULL;
    370   struct sln_tree* tree2 = NULL;
    371 
    372   const char* filename = "tree.sln";
    373   FILE* fp = NULL;
    374 
    375   CHK(sln_tree_create(sln, tree_args, &tree1) == RES_OK);
    376 
    377   CHK(fp = fopen(filename, "w+"));
    378 
    379   wargs.file = fp;
    380   CHK(sln_tree_write(NULL, &wargs) == RES_BAD_ARG);
    381   CHK(sln_tree_write(tree1, NULL) == RES_BAD_ARG);
    382   wargs.file = NULL;
    383   CHK(sln_tree_write(tree1, &wargs) == RES_BAD_ARG);
    384   wargs.file = fp;
    385   CHK(sln_tree_write(tree1, &wargs) == RES_OK);
    386   rewind(fp);
    387 
    388   rargs.metadata = tree_args->metadata;
    389   rargs.lines = tree_args->lines;
    390   rargs.file = fp;
    391   CHK(sln_tree_read(NULL, &rargs, &tree2) == RES_BAD_ARG);
    392   CHK(sln_tree_read(sln, NULL, &tree2) == RES_BAD_ARG);
    393   rargs.metadata = NULL;
    394   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_BAD_ARG);
    395   rargs.metadata = tree_args->metadata;
    396   rargs.lines = NULL;
    397   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_BAD_ARG);
    398   rargs.lines = tree_args->lines;
    399   rargs.file = NULL;
    400   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_BAD_ARG);
    401   rargs.file = fp;
    402   CHK(sln_tree_read(sln, &rargs, NULL) == RES_BAD_ARG);
    403   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_OK);
    404   fclose(fp);
    405 
    406   check_tree_equality(tree1, tree2);
    407   CHK(sln_tree_ref_put(tree2) == RES_OK);
    408 
    409   wargs.file = NULL;
    410   wargs.filename = filename;
    411   CHK(sln_tree_write(tree1, &wargs) == RES_OK);
    412 
    413   rargs.file = NULL;
    414   rargs.filename = "nop";
    415   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_IO_ERR);
    416   rargs.filename = filename;
    417   CHK(sln_tree_read(sln, &rargs, &tree2) == RES_OK);
    418 
    419   check_tree_equality(tree1, tree2);
    420 
    421   CHK(sln_tree_ref_put(tree2) == RES_OK);
    422   CHK(sln_tree_ref_put(tree1) == RES_OK);
    423 }
    424 
    425 /*******************************************************************************
    426  * Test function
    427  ******************************************************************************/
    428 int
    429 main(int argc, char** argv)
    430 {
    431   struct sln_device_create_args dev_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
    432   struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT;
    433   struct sln_device* sln = NULL;
    434 
    435   struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
    436   struct shtr* shtr = NULL;
    437   struct shtr_line_list_load_args line_load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
    438   struct shtr_line_list* line_list = NULL;
    439   struct shtr_isotope_metadata* metadata = NULL;
    440 
    441   FILE* fp_lines = NULL;
    442   FILE* fp_mdata = NULL;
    443   (void)argc, (void)argv;
    444 
    445   /* Generate the file of the isotope metadata */
    446   CHK(fp_mdata = tmpfile());
    447   fprintf(fp_mdata, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n");
    448   write_shtr_molecule(fp_mdata, &g_H2O);
    449   write_shtr_molecule(fp_mdata, &g_CO2);
    450   write_shtr_molecule(fp_mdata, &g_O3);
    451   rewind(fp_mdata);
    452 
    453   /* Generate the file of lines */
    454   CHK(fp_lines = tmpfile());
    455   write_shtr_lines(fp_lines, g_lines, g_nlines);
    456   rewind(fp_lines);
    457 
    458   /* Load the isotope metadata and the lines */
    459   shtr_args.verbose = 1;
    460   CHK(shtr_create(&shtr_args, &shtr) == RES_OK);
    461   CHK(shtr_isotope_metadata_load_stream(shtr, fp_mdata, NULL, &metadata) == RES_OK);
    462 
    463   line_load_args.filename = "stream";
    464   line_load_args.file = fp_lines;
    465   CHK(shtr_line_list_load(shtr, &line_load_args, &line_list) == RES_OK);
    466 
    467   CHK(fclose(fp_lines) == 0);
    468   CHK(fclose(fp_mdata) == 0);
    469 
    470   dev_args.verbose = 1;
    471   CHK(sln_device_create(&dev_args, &sln) == RES_OK);
    472 
    473   /* Create the mixture */
    474   tree_args.metadata = metadata;
    475   tree_args.lines = line_list;
    476   tree_args.molecules[SHTR_H2O].concentration = 1.0/3.0;
    477   tree_args.molecules[SHTR_H2O].cutoff = 25;
    478   tree_args.molecules[SHTR_CO2].concentration = 1.0/3.0;
    479   tree_args.molecules[SHTR_CO2].cutoff = 50;
    480   tree_args.molecules[SHTR_O3 ].concentration = 1.0/3.0;
    481   tree_args.molecules[SHTR_O3 ].cutoff = 25;
    482   tree_args.pressure = 1;
    483   tree_args.temperature = 296;
    484 
    485   test_tree(sln, &tree_args, line_list);
    486   test_tree_serialization(sln, &tree_args);
    487 
    488   /* Test a tree with a non default arity */
    489   tree_args.arity = 13;
    490   test_tree(sln, &tree_args, line_list);
    491   test_tree_serialization(sln, &tree_args);
    492 
    493   CHK(sln_device_ref_put(sln) == RES_OK);
    494   CHK(shtr_ref_put(shtr) == RES_OK);
    495   CHK(shtr_line_list_ref_put(line_list) == RES_OK);
    496   CHK(shtr_isotope_metadata_ref_put(metadata) == RES_OK);
    497   CHK(mem_allocated_size() == 0);
    498   return 0;
    499 }