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