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