star-line

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

sln_get.c (11218B)


      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 #define _POSIX_C_SOURCE 200112L /* getopt */
     20 
     21 #include "sln.h"
     22 
     23 #include <rsys/cstr.h>
     24 #include <rsys/math.h>
     25 #include <rsys/mem_allocator.h>
     26 
     27 #include <unistd.h>
     28 
     29 /* The maximum depth of a tree is 63, which is greater than the number of
     30  * lines that could actually be managed. Indeed, with such a depth, assuming
     31  * that there is a maximum of one line per leaf, it would be possible to
     32  * structure up to 2^63 lines. The problem would then be memory usage well
     33  * before this maximum depth.
     34  *
     35  * In fact, this limit is dictated by the number of bits used to encode the
     36  * descent masks, i.e., on a 64-bit integer. */
     37 #define MAX_DEPTH 64
     38 
     39 enum child { LEFT, RIGHT };
     40 
     41 enum output_type {
     42   OUTPUT_TREE_DESCRIPTOR,
     43   OUTPUT_NODE_DESCRIPTOR,
     44   OUTPUT_NODE_MESH,
     45   OUTPUT_NODE_VALUE,
     46   OUTPUT_COUNT__
     47 };
     48 
     49 struct args {
     50   const char* tree; /* NULL <=> read from standard input */
     51 
     52   enum output_type output_type;
     53 
     54   double wavenumber; /* Wave number at which the spectrum is evaluated */
     55 
     56   /* Step of the tree descent.
     57    *
     58    * A bit set in one of the masks determines which side to descend to move from
     59    * the current depth (which corresponds to the number of the set bit) to the
     60    * next one. Since the tree is binary, only one of the two masks is needed, as
     61    * the other can be deduced from the selected mask and the current depth.
     62    * However, calculating both descent masks allows bugs in the analysis of
     63    * arguments to be detected. */
     64   uint64_t descent_lmask; /* Left mask */
     65   uint64_t descent_rmask; /* Right mask */
     66   unsigned depth; /* Current depth in the tree */
     67 
     68   /* Miscellaneous */
     69   int quit;
     70   int verbose;
     71 };
     72 #define ARGS_DEFAULT__ {NULL, OUTPUT_TREE_DESCRIPTOR, 0,0,0,0,0,0}
     73 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
     74 
     75 struct cmd {
     76   struct args args;
     77 
     78   struct sln_device* sln;
     79   struct sln_tree* tree;
     80 
     81   struct shtr* shtr;
     82 };
     83 #define CMD_NULL__ {0}
     84 static const struct cmd CMD_NULL = CMD_NULL__;
     85 
     86 /*******************************************************************************
     87  * Helper functions
     88  ******************************************************************************/
     89 static void
     90 usage(FILE* stream)
     91 {
     92   fprintf(stream,
     93 "usage: sln-get [-hlmnrv] [-L nlevels] [-R nlevels] [-w wavenumber] [tree]\n");
     94 }
     95 
     96 static void
     97 tree_descent
     98   (struct args* args,
     99    const enum child child, /* side of the tree descent */
    100    const unsigned nlevels) /* Number of levels to descend in the tree */
    101 {
    102   uint64_t mask = 0;
    103   ASSERT(args);
    104 
    105   if(nlevels == 0 || args->depth >= MAX_DEPTH) return; /* Nothing to descent */
    106 
    107   if(nlevels > 64) {
    108     /* The mask can encode up to 64 levels of descent. Thus, if the number of
    109      * levels is greater than or equal to 64, all bits of the mask are
    110      * activated. */
    111     mask = ~0lu;
    112   } else {
    113     /* Set mask bits from bit 0 to the nlevels^th bit */
    114     mask = BIT_U64(nlevels)-1;
    115   }
    116 
    117   /* Move the mask to the current depth */
    118   mask <<= args->depth;
    119 
    120   /* Here we go: descent the tree */
    121   switch(child) {
    122     case LEFT:  args->descent_lmask |= mask; break;
    123     case RIGHT: args->descent_rmask |= mask; break;
    124     default: FATAL("Unreachable code\n"); break;
    125   }
    126   args->depth = MMIN(args->depth + nlevels, MAX_DEPTH);
    127 }
    128 
    129 static res_T
    130 args_init(struct args* args, int argc, char** argv)
    131 {
    132   unsigned nlvls = 0;
    133   int opt = 0;
    134   res_T res = RES_OK;
    135 
    136   ASSERT(args);
    137 
    138   *args = ARGS_DEFAULT;
    139 
    140   while((opt = getopt(argc, argv, "hL:lmnR:rvw:")) != -1) {
    141     switch(opt) {
    142       case 'h':
    143         usage(stdout);
    144         args->quit = 1;
    145         goto exit;
    146       case 'L':
    147         if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) {
    148           tree_descent(args, LEFT, nlvls);
    149         }
    150         break;
    151       case 'l': tree_descent(args, LEFT, 1); break;
    152       case 'R':
    153         if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) {
    154           tree_descent(args, RIGHT, nlvls);
    155         }
    156         break;
    157       case 'r': tree_descent(args, RIGHT, 1); break;
    158       case 'm': args->output_type = OUTPUT_NODE_MESH; break;
    159       case 'n': args->output_type = OUTPUT_NODE_DESCRIPTOR; break;
    160       case 'v': args->verbose += (args->verbose < 3); break;
    161       case 'w':
    162         args->output_type = OUTPUT_NODE_VALUE;
    163         res = cstr_to_double(optarg, &args->wavenumber);
    164         break;
    165       default: res = RES_BAD_ARG; break;
    166     }
    167     if(res != RES_OK) {
    168       if(optarg) {
    169         fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
    170           argv[0], optarg, opt);
    171       }
    172       goto error;
    173     }
    174   }
    175 
    176 #ifndef NDEBUG
    177   {
    178     const uint64_t descent_mask =
    179       args->depth >= 64 ? ~0ul : BIT_U64(args->depth)-1;
    180     ASSERT((args->descent_lmask | args->descent_rmask) == descent_mask);
    181   }
    182 #endif
    183 
    184   if(optind < argc) args->tree = argv[optind];
    185 
    186 exit:
    187   return res;
    188 error:
    189   usage(stderr);
    190   goto exit;
    191 }
    192 
    193 static void
    194 cmd_release(struct cmd* cmd)
    195 {
    196   ASSERT(cmd);
    197   if(cmd->sln) SLN(device_ref_put(cmd->sln));
    198   if(cmd->tree) SLN(tree_ref_put(cmd->tree));
    199   if(cmd->shtr) SHTR(ref_put(cmd->shtr));
    200 }
    201 
    202 static res_T
    203 cmd_init(struct cmd* cmd, const struct args* args)
    204 {
    205   struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
    206   struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL;
    207   struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
    208   res_T res = RES_OK;
    209 
    210   ASSERT(cmd && args);
    211 
    212   *cmd = CMD_NULL;
    213 
    214   shtr_args.verbose = args->verbose;
    215   res = shtr_create(&shtr_args, &cmd->shtr);
    216   if(res != RES_OK) goto error;
    217 
    218   sln_args.verbose = args->verbose;
    219   res = sln_device_create(&sln_args, &cmd->sln);
    220   if(res != RES_OK) goto error;
    221 
    222   tree_args.shtr = cmd->shtr;
    223   if(args->tree) {
    224     tree_args.file = NULL;
    225     tree_args.filename = args->tree;
    226   } else {
    227     tree_args.file = stdin;
    228     tree_args.filename = "stdin";
    229   }
    230   res = sln_tree_read(cmd->sln, &tree_args, &cmd->tree);
    231   if(res != RES_OK) goto error;
    232 
    233   cmd->args = *args;
    234 
    235 exit:
    236   return res;
    237 error:
    238   cmd_release(cmd);
    239   *cmd = CMD_NULL;
    240   goto exit;
    241 }
    242 
    243 static res_T
    244 print_tree_descriptor(const struct cmd* cmd)
    245 {
    246   struct sln_tree_desc desc = SLN_TREE_DESC_NULL;
    247   res_T res = RES_OK;
    248   ASSERT(cmd);
    249 
    250   res = sln_tree_get_desc(cmd->tree, &desc);
    251   if(res != RES_OK) goto error;
    252 
    253   printf("#lines:           %lu\n", (unsigned long)desc.nlines);
    254   printf("#nodes:           %lu\n", (unsigned long)desc.nnodes);
    255   printf("tree depth:       %u\n", desc.depth);
    256   printf("#vertices:        %lu\n", (unsigned long)desc.nvertices);
    257   printf("type:             %s\n", sln_mesh_type_cstr(desc.mesh_type));
    258   printf("decimation error: %.4e\n", desc.mesh_decimation_err);
    259   printf("line profile:     %s\n", sln_line_profile_cstr(desc.line_profile));
    260   printf("#lines per leaf:  %lu\n", (unsigned long)desc.max_nlines_per_leaf);
    261 
    262 exit:
    263   return res;
    264 error:
    265   goto exit;
    266 }
    267 
    268 static const struct sln_node* /* NULL <=> tree is empty */
    269 get_node(const struct cmd* cmd, unsigned* node_depth/*can be NULL*/)
    270 {
    271   const struct sln_node* node = NULL;
    272   unsigned depth = 0;
    273   unsigned i = 0;
    274   ASSERT(cmd);
    275 
    276   node = sln_tree_get_root(cmd->tree);
    277   if(node == NULL) goto exit; /* Tree is empty */
    278 
    279   FOR_EACH(i, 0, cmd->args.depth) {
    280     const uint64_t mask = BIT_U64(i);
    281 
    282     if(sln_node_is_leaf(node)) break;
    283 
    284     if(mask & cmd->args.descent_lmask) {
    285       node = sln_node_get_child(node, 0);
    286     } else if(mask & cmd->args.descent_rmask) {
    287       node = sln_node_get_child(node, 1);
    288     } else {
    289       FATAL("Unreachable code\n");
    290     }
    291     ++depth;
    292   }
    293 
    294 exit:
    295   if(node_depth) *node_depth = depth;
    296   return node;
    297 }
    298 
    299 static res_T
    300 print_node_descriptor(const struct cmd* cmd)
    301 {
    302   const struct sln_node* node = NULL;
    303   struct sln_node_desc desc = SLN_NODE_DESC_NULL;
    304   unsigned depth = 0;
    305   res_T res = RES_OK;
    306   ASSERT(cmd);
    307 
    308   if((node = get_node(cmd, &depth)) == NULL) goto exit; /* tree is empty */
    309 
    310   res = sln_node_get_desc(node, &desc);
    311   if(res != RES_OK) goto error;
    312 
    313   printf("level:     %u\n", depth);
    314   printf("#lines:    %lu\n", (unsigned long)desc.nlines);
    315   printf("#vertices: %lu\n", (unsigned long)desc.nvertices);
    316 
    317 exit:
    318   return res;
    319 error:
    320   goto exit;
    321 }
    322 
    323 static res_T
    324 print_mesh(const struct cmd* cmd)
    325 {
    326   struct sln_mesh mesh = SLN_MESH_NULL;
    327   const struct sln_node* node = NULL;
    328   size_t i = 0;
    329   res_T res = RES_OK;
    330   ASSERT(cmd);
    331 
    332   if((node = get_node(cmd, NULL)) == NULL) goto exit; /* tree is empty */
    333 
    334   res = sln_node_get_mesh(cmd->tree, node, &mesh);
    335   if(res != RES_OK) goto error;
    336 
    337   FOR_EACH(i, 0, mesh.nvertices) {
    338     printf("%g %g\n",
    339       mesh.vertices[i].wavenumber,
    340       mesh.vertices[i].ka);
    341   }
    342 
    343 exit:
    344   return res;
    345 error:
    346   goto exit;
    347 }
    348 
    349 static res_T
    350 print_node_value(const struct cmd* cmd)
    351 {
    352   struct sln_mesh mesh = SLN_MESH_NULL;
    353   const struct sln_node* node = NULL;
    354   double val_mesh = 0;
    355   double val_node = 0;
    356   res_T res = RES_OK;
    357   ASSERT(cmd);
    358 
    359   if((node = get_node(cmd, NULL)) == NULL) goto exit; /* tree is empty */
    360 
    361   res = sln_node_get_mesh(cmd->tree, node, &mesh);
    362   if(res != RES_OK) goto error;
    363 
    364   val_mesh = sln_mesh_eval(&mesh, cmd->args.wavenumber);
    365   val_node = sln_node_eval(cmd->tree, node, cmd->args.wavenumber);
    366 
    367   printf("ka(%e) = %e ~ %e\n",  cmd->args.wavenumber, val_node, val_mesh);
    368 
    369 exit:
    370   return res;
    371 error:
    372   goto exit;
    373 }
    374 
    375 static res_T
    376 cmd_run(const struct cmd* cmd)
    377 {
    378   res_T res = RES_OK;
    379 
    380   switch(cmd->args.output_type) {
    381     case OUTPUT_TREE_DESCRIPTOR:
    382       res = print_tree_descriptor(cmd);
    383       break;
    384     case OUTPUT_NODE_DESCRIPTOR:
    385       res = print_node_descriptor(cmd);
    386       break;
    387     case OUTPUT_NODE_MESH:
    388       res = print_mesh(cmd);
    389       break;
    390     case OUTPUT_NODE_VALUE:
    391       res = print_node_value(cmd);
    392       break;
    393     default: FATAL("Unreachable code\n"); break;
    394   }
    395   if(res != RES_OK) goto error;
    396 
    397 exit:
    398   return res;
    399 error:
    400   goto exit;
    401 }
    402 
    403 /*******************************************************************************
    404  * The program
    405  ******************************************************************************/
    406 int
    407 main(int argc, char** argv)
    408 {
    409   struct args args = ARGS_DEFAULT;
    410   struct cmd cmd = CMD_NULL;
    411   int err = 0;
    412   res_T res = RES_OK;
    413 
    414   if((res = args_init(&args, argc, argv)) != RES_OK) goto error;
    415   if(args.quit) goto exit;
    416 
    417   if((res = cmd_init(&cmd, &args)) != RES_OK) goto error;
    418   if((res = cmd_run(&cmd)) != RES_OK) goto error;
    419 
    420 exit:
    421   cmd_release(&cmd);
    422   CHK(mem_allocated_size() == 0);
    423   return err;
    424 error:
    425   err = 1;
    426   goto exit;
    427 }