sln_get.c (14946B)
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_LEVEL_DESCRIPTOR, 43 OUTPUT_NODE_DESCRIPTOR, 44 OUTPUT_NODE_MESH, 45 OUTPUT_NODE_VALUE, 46 OUTPUT_TREE_DESCRIPTOR, 47 OUTPUT_COUNT__ 48 }; 49 50 struct args { 51 const char* tree; /* NULL <=> read from standard input */ 52 const char* molparams; 53 const char* lines; 54 55 enum output_type output_type; 56 57 double wavenumber; /* Wave number at which the spectrum is evaluated */ 58 59 /* Step of the tree descent. 60 * 61 * A bit set in one of the masks determines which side to descend to move from 62 * the current depth (which corresponds to the number of the set bit) to the 63 * next one. Since the tree is binary, only one of the two masks is needed, as 64 * the other can be deduced from the selected mask and the current depth. 65 * However, calculating both descent masks allows bugs in the analysis of 66 * arguments to be detected. */ 67 uint64_t descent_lmask; /* Left mask */ 68 uint64_t descent_rmask; /* Right mask */ 69 unsigned depth; /* Current depth in the tree */ 70 71 /* Miscellaneous */ 72 unsigned level; /* Queried level */ 73 int quit; 74 int verbose; 75 int lines_in_shtr_format; 76 }; 77 #define ARGS_DEFAULT__ {NULL,NULL,NULL,OUTPUT_TREE_DESCRIPTOR,0,0,0,0,0,0,0,0} 78 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 79 80 struct cmd { 81 struct args args; 82 83 struct sln_device* sln; 84 struct sln_tree* tree; 85 86 struct shtr* shtr; 87 struct shtr_isotope_metadata* molparams; 88 struct shtr_line_list* lines; 89 }; 90 #define CMD_NULL__ {0} 91 static const struct cmd CMD_NULL = CMD_NULL__; 92 93 /******************************************************************************* 94 * Helper functions 95 ******************************************************************************/ 96 static void 97 usage(FILE* stream) 98 { 99 fprintf(stream, 100 "usage: sln-get [-hlmnrsv] [-d level] [-L nlevels] [-R nlevels] [-w wavenumber]\n" 101 " -i lines -p molparams [tree]\n"); 102 } 103 104 static void 105 tree_descent 106 (struct args* args, 107 const enum child child, /* side of the tree descent */ 108 const unsigned nlevels) /* Number of levels to descend in the tree */ 109 { 110 uint64_t mask = 0; 111 ASSERT(args); 112 113 if(nlevels == 0 || args->depth >= MAX_DEPTH) return; /* Nothing to descent */ 114 115 if(nlevels > 64) { 116 /* The mask can encode up to 64 levels of descent. Thus, if the number of 117 * levels is greater than or equal to 64, all bits of the mask are 118 * activated. */ 119 mask = ~0lu; 120 } else { 121 /* Set mask bits from bit 0 to the nlevels^th bit */ 122 mask = BIT_U64(nlevels)-1; 123 } 124 125 /* Move the mask to the current depth */ 126 mask <<= args->depth; 127 128 /* Here we go: descent the tree */ 129 switch(child) { 130 case LEFT: args->descent_lmask |= mask; break; 131 case RIGHT: args->descent_rmask |= mask; break; 132 default: FATAL("Unreachable code\n"); break; 133 } 134 args->depth = MMIN(args->depth + nlevels, MAX_DEPTH); 135 } 136 137 static res_T 138 args_init(struct args* args, int argc, char** argv) 139 { 140 unsigned nlvls = 0; 141 int opt = 0; 142 res_T res = RES_OK; 143 144 ASSERT(args); 145 146 *args = ARGS_DEFAULT; 147 148 while((opt = getopt(argc, argv, "d:hi:L:lmnp:R:rsvw:")) != -1) { 149 switch(opt) { 150 case 'd': 151 args->output_type = OUTPUT_LEVEL_DESCRIPTOR; 152 res = cstr_to_uint(optarg, &args->level); 153 break; 154 case 'h': 155 usage(stdout); 156 args->quit = 1; 157 goto exit; 158 case 'i': args->lines = optarg; break; 159 case 'L': 160 if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) { 161 tree_descent(args, LEFT, nlvls); 162 } 163 break; 164 case 'l': tree_descent(args, LEFT, 1); break; 165 case 'm': args->output_type = OUTPUT_NODE_MESH; break; 166 case 'n': args->output_type = OUTPUT_NODE_DESCRIPTOR; break; 167 case 'p': args->molparams = optarg; break; 168 case 'R': 169 if((res = cstr_to_uint(optarg, &nlvls)) == RES_OK) { 170 tree_descent(args, RIGHT, nlvls); 171 } 172 break; 173 case 'r': tree_descent(args, RIGHT, 1); break; 174 case 's': args->lines_in_shtr_format = 1; break; 175 case 'v': args->verbose += (args->verbose < 3); break; 176 case 'w': 177 args->output_type = OUTPUT_NODE_VALUE; 178 res = cstr_to_double(optarg, &args->wavenumber); 179 break; 180 default: res = RES_BAD_ARG; break; 181 } 182 if(res != RES_OK) { 183 if(optarg) { 184 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 185 argv[0], optarg, opt); 186 } 187 goto error; 188 } 189 } 190 191 #define MANDATORY(Cond, Name, Opt) { \ 192 if(!(Cond)) { \ 193 fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \ 194 res = RES_BAD_ARG; \ 195 goto error; \ 196 } \ 197 } (void)0 198 MANDATORY(args->molparams, "molparams", 'p'); 199 MANDATORY(args->lines, "line list", 'i'); 200 #undef MANDATORY 201 202 #ifndef NDEBUG 203 { 204 const uint64_t descent_mask = 205 args->depth >= 64 ? ~0ul : BIT_U64(args->depth)-1; 206 ASSERT((args->descent_lmask | args->descent_rmask) == descent_mask); 207 } 208 #endif 209 210 if(optind < argc) args->tree = argv[optind]; 211 212 exit: 213 return res; 214 error: 215 usage(stderr); 216 goto exit; 217 } 218 219 static void 220 cmd_release(struct cmd* cmd) 221 { 222 ASSERT(cmd); 223 if(cmd->sln) SLN(device_ref_put(cmd->sln)); 224 if(cmd->tree) SLN(tree_ref_put(cmd->tree)); 225 if(cmd->shtr) SHTR(ref_put(cmd->shtr)); 226 if(cmd->molparams) SHTR(isotope_metadata_ref_put(cmd->molparams)); 227 if(cmd->lines) SHTR(line_list_ref_put(cmd->lines)); 228 } 229 230 static res_T 231 load_lines(struct cmd* cmd, const struct args* args) 232 { 233 res_T res = RES_OK; 234 ASSERT(cmd && args); 235 236 if(args->lines_in_shtr_format) { 237 struct shtr_line_list_read_args read_args = SHTR_LINE_LIST_READ_ARGS_NULL; 238 239 /* Loads lines from data serialized by the Star-HITRAN library */ 240 read_args.filename = args->lines; 241 res = shtr_line_list_read(cmd->shtr, &read_args, &cmd->lines); 242 if(res != RES_OK) goto error; 243 244 } else { 245 struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL; 246 247 /* Loads lines from a file in HITRAN format */ 248 load_args.filename = args->lines; 249 res = shtr_line_list_load(cmd->shtr, &load_args, &cmd->lines); 250 if(res != RES_OK) goto error; 251 } 252 253 exit: 254 return res; 255 error: 256 if(cmd->lines) { SHTR(line_list_ref_put(cmd->lines)); cmd->lines = NULL; } 257 goto exit; 258 } 259 260 static res_T 261 cmd_init(struct cmd* cmd, const struct args* args) 262 { 263 struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 264 struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL; 265 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 266 res_T res = RES_OK; 267 268 ASSERT(cmd && args); 269 270 *cmd = CMD_NULL; 271 272 shtr_args.verbose = args->verbose; 273 res = shtr_create(&shtr_args, &cmd->shtr); 274 if(res != RES_OK) goto error; 275 276 res = shtr_isotope_metadata_load(cmd->shtr, args->molparams, &cmd->molparams); 277 if(res != RES_OK) goto error; 278 279 res = load_lines(cmd, args); 280 if(res != RES_OK) goto error; 281 282 sln_args.verbose = args->verbose; 283 res = sln_device_create(&sln_args, &cmd->sln); 284 if(res != RES_OK) goto error; 285 286 tree_args.metadata = cmd->molparams; 287 tree_args.lines = cmd->lines; 288 if(args->tree) { 289 tree_args.file = NULL; 290 tree_args.filename = args->tree; 291 } else { 292 tree_args.file = stdin; 293 tree_args.filename = "stdin"; 294 } 295 res = sln_tree_read(cmd->sln, &tree_args, &cmd->tree); 296 if(res != RES_OK) goto error; 297 298 cmd->args = *args; 299 300 exit: 301 return res; 302 error: 303 cmd_release(cmd); 304 *cmd = CMD_NULL; 305 goto exit; 306 } 307 308 static res_T 309 print_level_descriptor(const struct cmd* cmd) 310 { 311 /* Stack for visiting the tree depth-first */ 312 struct { const struct sln_node* node; unsigned level; } stack[MAX_DEPTH]; 313 int istack = 0; 314 315 /* Node data */ 316 struct sln_node_desc desc = SLN_NODE_DESC_NULL; 317 const struct sln_node* node = NULL; 318 319 /* Level descriptor */ 320 size_t nvertices = 0; 321 size_t nnodes = 0; 322 323 /* Miscellaneous */ 324 unsigned level = 0; 325 res_T res = RES_OK; 326 327 ASSERT(cmd); /* Precondition */ 328 329 /* Push a dummy node which, once pop up, whill mark the end of recursion */ 330 stack[istack].node = NULL; 331 stack[istack].level = UINT_MAX; 332 ++istack; 333 334 node = sln_tree_get_root(cmd->tree); 335 336 while(node) { 337 ASSERT(level <= cmd->args.level); 338 339 if(!sln_node_is_leaf(node) && level < cmd->args.level) { 340 /* Continue down the tree */ 341 ++level; 342 343 /* Push the right child */ 344 stack[istack ].node = sln_node_get_child(node, 1); 345 stack[istack++].level = level; 346 347 node = sln_node_get_child(node, 0); /* Visit the left child */ 348 349 } else { 350 /* The queried level or a leaf is reached, update the descriptor */ 351 if((res = sln_node_get_desc(node, &desc)) != RES_OK) goto error; 352 nvertices += desc.nvertices; 353 ++nnodes; 354 355 /* Pop the next node */; 356 node = stack[--istack].node; 357 level = stack[ istack].level; 358 } 359 } 360 361 /* Print the level description */ 362 printf("#nodes: %lu\n", (unsigned long)nnodes); 363 printf("#vertices: %lu\n", (unsigned long)nvertices); 364 365 exit: 366 return res; 367 error: 368 goto exit; 369 } 370 371 static const struct sln_node* /* NULL <=> tree is empty */ 372 get_node(const struct cmd* cmd, unsigned* node_depth/*can be NULL*/) 373 { 374 const struct sln_node* node = NULL; 375 unsigned depth = 0; 376 unsigned i = 0; 377 ASSERT(cmd); 378 379 node = sln_tree_get_root(cmd->tree); 380 if(node == NULL) goto exit; /* Tree is empty */ 381 382 FOR_EACH(i, 0, cmd->args.depth) { 383 const uint64_t mask = BIT_U64(i); 384 385 if(sln_node_is_leaf(node)) break; 386 387 if(mask & cmd->args.descent_lmask) { 388 node = sln_node_get_child(node, 0); 389 } else if(mask & cmd->args.descent_rmask) { 390 node = sln_node_get_child(node, 1); 391 } else { 392 FATAL("Unreachable code\n"); 393 } 394 ++depth; 395 } 396 397 exit: 398 if(node_depth) *node_depth = depth; 399 return node; 400 } 401 402 static res_T 403 print_node_descriptor(const struct cmd* cmd) 404 { 405 const struct sln_node* node = NULL; 406 struct sln_node_desc desc = SLN_NODE_DESC_NULL; 407 unsigned depth = 0; 408 res_T res = RES_OK; 409 ASSERT(cmd); 410 411 if((node = get_node(cmd, &depth)) == NULL) goto exit; /* tree is empty */ 412 413 res = sln_node_get_desc(node, &desc); 414 if(res != RES_OK) goto error; 415 416 printf("level: %u\n", depth); 417 printf("#lines: %lu\n", (unsigned long)desc.nlines); 418 printf("#vertices: %lu\n", (unsigned long)desc.nvertices); 419 420 exit: 421 return res; 422 error: 423 goto exit; 424 } 425 426 static res_T 427 print_mesh(const struct cmd* cmd) 428 { 429 struct sln_mesh mesh = SLN_MESH_NULL; 430 const struct sln_node* node = NULL; 431 size_t i = 0; 432 res_T res = RES_OK; 433 ASSERT(cmd); 434 435 if((node = get_node(cmd, NULL)) == NULL) goto exit; /* tree is empty */ 436 437 res = sln_node_get_mesh(cmd->tree, node, &mesh); 438 if(res != RES_OK) goto error; 439 440 FOR_EACH(i, 0, mesh.nvertices) { 441 printf("%g %g\n", 442 mesh.vertices[i].wavenumber, 443 mesh.vertices[i].ka); 444 } 445 446 exit: 447 return res; 448 error: 449 goto exit; 450 } 451 452 static res_T 453 print_node_value(const struct cmd* cmd) 454 { 455 struct sln_mesh mesh = SLN_MESH_NULL; 456 const struct sln_node* node = NULL; 457 double val_mesh = 0; 458 double val_node = 0; 459 res_T res = RES_OK; 460 ASSERT(cmd); 461 462 if((node = get_node(cmd, NULL)) == NULL) goto exit; /* tree is empty */ 463 464 res = sln_node_get_mesh(cmd->tree, node, &mesh); 465 if(res != RES_OK) goto error; 466 467 val_mesh = sln_mesh_eval(&mesh, cmd->args.wavenumber); 468 val_node = sln_node_eval(cmd->tree, node, cmd->args.wavenumber); 469 470 printf("ka(%e) = %e ~ %e\n", cmd->args.wavenumber, val_node, val_mesh); 471 472 exit: 473 return res; 474 error: 475 goto exit; 476 } 477 478 static res_T 479 print_tree_descriptor(const struct cmd* cmd) 480 { 481 struct sln_tree_desc desc = SLN_TREE_DESC_NULL; 482 res_T res = RES_OK; 483 ASSERT(cmd); 484 485 res = sln_tree_get_desc(cmd->tree, &desc); 486 if(res != RES_OK) goto error; 487 488 printf("#lines: %lu\n", (unsigned long)desc.nlines); 489 printf("#nodes: %lu\n", (unsigned long)desc.nnodes); 490 printf("tree depth: %u\n", desc.depth); 491 printf("#vertices: %lu\n", (unsigned long)desc.nvertices); 492 printf("type: %s\n", sln_mesh_type_cstr(desc.mesh_type)); 493 printf("decimation error: %.4e\n", desc.mesh_decimation_err); 494 printf("line profile: %s\n", sln_line_profile_cstr(desc.line_profile)); 495 printf("#lines per leaf: %lu\n", (unsigned long)desc.max_nlines_per_leaf); 496 497 exit: 498 return res; 499 error: 500 goto exit; 501 } 502 503 static res_T 504 cmd_run(const struct cmd* cmd) 505 { 506 res_T res = RES_OK; 507 508 switch(cmd->args.output_type) { 509 case OUTPUT_LEVEL_DESCRIPTOR: 510 res = print_level_descriptor(cmd); 511 break; 512 case OUTPUT_NODE_DESCRIPTOR: 513 res = print_node_descriptor(cmd); 514 break; 515 case OUTPUT_NODE_MESH: 516 res = print_mesh(cmd); 517 break; 518 case OUTPUT_NODE_VALUE: 519 res = print_node_value(cmd); 520 break; 521 case OUTPUT_TREE_DESCRIPTOR: 522 res = print_tree_descriptor(cmd); 523 break; 524 default: FATAL("Unreachable code\n"); break; 525 } 526 if(res != RES_OK) goto error; 527 528 exit: 529 return res; 530 error: 531 goto exit; 532 } 533 534 /******************************************************************************* 535 * The program 536 ******************************************************************************/ 537 int 538 main(int argc, char** argv) 539 { 540 struct args args = ARGS_DEFAULT; 541 struct cmd cmd = CMD_NULL; 542 int err = 0; 543 res_T res = RES_OK; 544 545 if((res = args_init(&args, argc, argv)) != RES_OK) goto error; 546 if(args.quit) goto exit; 547 548 if((res = cmd_init(&cmd, &args)) != RES_OK) goto error; 549 if((res = cmd_run(&cmd)) != RES_OK) goto error; 550 551 exit: 552 cmd_release(&cmd); 553 CHK(mem_allocated_size() == 0); 554 return err; 555 error: 556 err = 1; 557 goto exit; 558 }