sln_tree.c (27572B)
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 "sln.h" 20 #include "sln_device_c.h" 21 #include "sln_line.h" 22 #include "sln_tree_c.h" 23 24 #include <star/shtr.h> 25 #include <star/ssp.h> 26 27 #include <rsys/algorithm.h> 28 #include <rsys/cstr.h> 29 #include <rsys/math.h> 30 31 #include <omp.h> 32 33 struct stream { 34 const char* name; 35 FILE* fp; 36 int intern_fp; /* Define if the stream was internally opened */ 37 }; 38 static const struct stream STREAM_NULL = {NULL, NULL, 0}; 39 40 /******************************************************************************* 41 * Helper functions 42 ******************************************************************************/ 43 /* Check that the sum of the molecular concentrations is equal to 1 */ 44 static res_T 45 check_molecule_concentration 46 (struct sln_device* sln, 47 const char* caller, 48 const struct sln_tree_create_args* args) 49 { 50 int i = 0; 51 double sum = 0; 52 ASSERT(sln && caller && args); 53 54 FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) { 55 if(i == SHTR_MOLECULE_ID_NULL) continue; 56 sum += args->molecules[i].concentration; 57 } 58 59 /* The sum of molecular concentrations must be less than or equal to 1. It may 60 * be less than 1 if the remaining part of the mixture is (implicitly) defined 61 * as a radiatively inactive gas */ 62 if(sum > 1 && sum-1 > 1e-6) { 63 ERROR(sln, 64 "%s: the sum of molecule concentrations is greater than 1: %g\n", 65 caller, sum); 66 return RES_BAD_ARG; 67 } 68 69 return RES_OK; 70 } 71 72 /* Verify that the isotope abundance are valids */ 73 static res_T 74 check_molecule_isotope_abundances 75 (struct sln_device* sln, 76 const char* caller, 77 const struct sln_molecule* molecule) 78 { 79 int i = 0; 80 double sum = 0; 81 ASSERT(sln && caller && molecule); 82 83 /* The isotopic abundances are the default ones. Nothing to do */ 84 if(!molecule->non_default_isotope_abundances) return RES_OK; 85 86 /* The isotopic abundances are not the default ones. 87 * Verify that they are valid ... */ 88 FOR_EACH(i, 0, SHTR_MAX_ISOTOPE_COUNT) { 89 if(molecule->isotopes[i].abundance < 0) { 90 const int isotope_id = i + 1; /* isotope id in [1, 10] */ 91 ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n", 92 caller, isotope_id, shtr_molecule_cstr(i), 93 molecule->isotopes[i].abundance); 94 return RES_BAD_ARG; 95 } 96 97 sum += molecule->isotopes[i].abundance; 98 } 99 100 /* ... and that their sum equals 1 */ 101 if(!eq_eps(sum, 1, 1e-6)) { 102 ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n", 103 caller, shtr_molecule_cstr(i), sum); 104 return RES_BAD_ARG; 105 } 106 107 return RES_OK; 108 } 109 110 static res_T 111 check_molecules 112 (struct sln_device* sln, 113 const char* caller, 114 const struct sln_tree_create_args* args) 115 { 116 char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0}; 117 118 double concentrations_sum = 0; 119 size_t iline = 0; 120 size_t nlines = 0; 121 res_T res = RES_OK; 122 ASSERT(args->lines); 123 124 res = check_molecule_concentration(sln, caller, args); 125 if(res != RES_OK) return res; 126 127 /* Iterate over the lines to define which molecules has to be checked, i.e., 128 * the ones used in the mixture */ 129 SHTR(line_list_get_size(args->lines, &nlines)); 130 FOR_EACH(iline, 0, nlines) { 131 struct shtr_line line = SHTR_LINE_NULL; 132 const struct sln_molecule* molecule = NULL; 133 134 SHTR(line_list_at(args->lines, iline, &line)); 135 136 /* This molecule was already checked */ 137 if(molecule_ok[line.molecule_id]) continue; 138 139 molecule = args->molecules + line.molecule_id; 140 141 if(molecule->concentration == 0) { 142 /* A molecular concentration of zero is allowed, but may be a user error, 143 * as 0 is the default concentration in the tree creation arguments. 144 * Therefore, warn the user about this value so that they can determine 145 * whether or not it is an error on their part. */ 146 WARN(sln, "%s: the concentration of %s is zero.\n", 147 caller, shtr_molecule_cstr(line.molecule_id)); 148 149 } else if(molecule->concentration < 0) { 150 /* Concentration cannot be negative... */ 151 ERROR(sln, "%s: invalid %s concentration: %g.\n", 152 FUNC_NAME, shtr_molecule_cstr(line.molecule_id), 153 molecule->concentration); 154 return RES_BAD_ARG; 155 } 156 157 concentrations_sum += molecule->concentration; 158 159 if(molecule->cutoff <= 0) { 160 /* ... cutoff either */ 161 ERROR(sln, "%s: invalid %s cutoff: %g.\n", 162 caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff); 163 return RES_BAD_ARG; 164 } 165 166 res = check_molecule_isotope_abundances(sln, caller, molecule); 167 if(res != RES_OK) return res; 168 169 molecule_ok[line.molecule_id] = 1; 170 } 171 172 /* The sum of molecular concentrations must be less than or equal to 1. It may 173 * be less than 1 if the remaining part of the mixture is (implicitly) defined 174 * as a radiatively inactive gas */ 175 if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) { 176 ERROR(sln, 177 "%s: the sum of molecule concentrations is greater than 1: %g\n", 178 caller, concentrations_sum); 179 return RES_BAD_ARG; 180 } 181 182 return RES_OK; 183 } 184 185 static res_T 186 check_sln_tree_create_args 187 (struct sln_device* sln, 188 const char* caller, 189 const struct sln_tree_create_args* args) 190 { 191 res_T res = RES_OK; 192 ASSERT(sln && caller); 193 194 if(!args) return RES_BAD_ARG; 195 196 if(!args->metadata) { 197 ERROR(sln, "%s: the isotope metadata are missing.\n", caller); 198 return RES_BAD_ARG; 199 } 200 201 if(!args->lines) { 202 ERROR(sln, "%s: the list of lines is missing.\n", caller); 203 return RES_BAD_ARG; 204 } 205 206 if(args->nvertices_hint == 0) { 207 ERROR(sln, 208 "%s: invalid hint on the number of vertices around the line center %lu.\n", 209 caller, (unsigned long)args->nvertices_hint); 210 return RES_BAD_ARG; 211 } 212 213 if(args->mesh_decimation_err < 0) { 214 ERROR(sln, "%s: invalid decimation error %g.\n", 215 caller, args->mesh_decimation_err); 216 return RES_BAD_ARG; 217 } 218 219 if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) { 220 ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type); 221 return RES_BAD_ARG; 222 } 223 224 if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) { 225 ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile); 226 return RES_BAD_ARG; 227 } 228 229 if(args->arity < 2 || args->arity > SLN_TREE_ARITY_MAX) { 230 ERROR(sln, "%s: invalid arity %u. It must be in [2, %d]\n", 231 caller, args->arity, SLN_TREE_ARITY_MAX); 232 return RES_BAD_ARG; 233 } 234 235 if(args->leaf_nlines < 1 || args->leaf_nlines > SLN_LEAF_NLINES_MAX) { 236 ERROR(sln, "%s: invalid number of lines per leaf %u. It must be in [1, %d]\n", 237 caller, args->leaf_nlines, SLN_LEAF_NLINES_MAX); 238 return RES_BAD_ARG; 239 } 240 241 if(args->nthreads_hint == 0) { 242 ERROR(sln, "%s: invalid number of threads %u\n", 243 caller, args->nthreads_hint); 244 return RES_BAD_ARG; 245 } 246 247 res = check_molecules(sln, caller, args); 248 if(res != RES_OK) return res; 249 250 return RES_OK; 251 } 252 253 static res_T 254 check_sln_tree_read_args 255 (struct sln_device* sln, 256 const char* caller, 257 const struct sln_tree_read_args* args) 258 { 259 if(!args) return RES_BAD_ARG; 260 261 if(!args->metadata) { 262 ERROR(sln, "%s: the isotope metadata are missing.\n", caller); 263 return RES_BAD_ARG; 264 } 265 266 if(!args->lines) { 267 ERROR(sln, "%s: the list of lines is missing.\n", caller); 268 return RES_BAD_ARG; 269 } 270 271 if(!args->file && !args->filename) { 272 ERROR(sln, 273 "%s: the source file is missing. No file name or stream is provided.\n", 274 caller); 275 return RES_BAD_ARG; 276 } 277 278 return RES_OK; 279 } 280 281 static res_T 282 check_sln_tree_write_args 283 (struct sln_device* sln, 284 const char* caller, 285 const struct sln_tree_write_args* args) 286 { 287 if(!args) return RES_BAD_ARG; 288 289 if(!args->file && !args->filename) { 290 ERROR(sln, 291 "%s: the destination file is missing. " 292 "No file name or stream is provided.\n", 293 caller); 294 return RES_BAD_ARG; 295 } 296 297 return RES_OK; 298 } 299 static INLINE void 300 stream_release(struct stream* stream) 301 { 302 ASSERT(stream); 303 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 304 } 305 306 static res_T 307 stream_init 308 (struct sln_device* sln, 309 const char* caller, 310 const char* name, /* NULL <=> default stream name */ 311 FILE* fp, /* NULL <=> open file "name" */ 312 const char* mode, /* mode in fopen */ 313 struct stream* stream) 314 { 315 res_T res = RES_OK; 316 317 ASSERT(sln && caller && stream); 318 ASSERT(fp || (name && mode)); 319 320 *stream = STREAM_NULL; 321 322 if(fp) { 323 stream->intern_fp = 0; 324 stream->name = name ? name : "stream"; 325 stream->fp = fp; 326 327 } else { 328 stream->intern_fp = 1; 329 stream->name = name; 330 if(!(stream->fp = fopen(name, mode))) { 331 ERROR(sln, "%s:%s: error opening file -- %s\n", 332 caller, name, strerror(errno)); 333 res = RES_IO_ERR; 334 goto error; 335 } 336 } 337 338 exit: 339 return res; 340 error: 341 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 342 goto exit; 343 } 344 345 static res_T 346 create_tree 347 (struct sln_device* sln, 348 const char* caller, 349 struct sln_tree** out_tree) 350 { 351 struct sln_tree* tree = NULL; 352 res_T res = RES_OK; 353 ASSERT(sln && caller && out_tree); 354 355 tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree)); 356 if(!tree) { 357 ERROR(sln, "%s: could not allocate the tree data structure.\n", 358 caller); 359 res = RES_MEM_ERR; 360 goto error; 361 } 362 ref_init(&tree->ref); 363 SLN(device_ref_get(sln)); 364 tree->sln = sln; 365 darray_node_init(sln->allocator, &tree->nodes); 366 darray_vertex_init(sln->allocator, &tree->vertices); 367 368 exit: 369 *out_tree = tree; 370 return res; 371 error: 372 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 373 goto exit; 374 } 375 376 static INLINE int 377 cmp_nu_vtx(const void* key, const void* item) 378 { 379 const float nu = *((const float*)key); 380 const struct sln_vertex* vtx = item; 381 if(nu < vtx->wavenumber) return -1; 382 if(nu > vtx->wavenumber) return +1; 383 return 0; 384 } 385 386 static res_T 387 build_node_cumulative 388 (const struct sln_tree* tree, 389 const struct sln_node* node, 390 const double nu, /* [cm^-1] */ 391 double proba[SLN_TREE_ARITY_MAX], 392 double cumul[SLN_TREE_ARITY_MAX]) 393 { 394 struct sln_mesh mesh = SLN_MESH_NULL; 395 double ka = 0; 396 double sum = 0; 397 unsigned i=0, n=0; 398 res_T res = RES_OK; 399 ASSERT(tree && node && proba && cumul); 400 401 n = sln_node_get_child_count(tree, node); 402 ASSERT(n <= SLN_TREE_ARITY_MAX); 403 404 FOR_EACH(i, 0, n) { 405 const struct sln_node* child = sln_node_get_child(tree, node, i); 406 407 SLN(node_get_mesh(tree, child, &mesh)); 408 ka = sln_mesh_eval(&mesh, nu); 409 410 sum += ka; 411 cumul[i] = sum; 412 proba[i] = ka; 413 } 414 415 /* No lines could be sampled because none of them affect the absorption at the 416 * given wave number */ 417 if(sum == 0) { 418 res = RES_BAD_ARG; 419 goto error; 420 } 421 422 /* Check the criterion of transition importance sampling, i.e. the value of 423 * the parent node must be greater than or equal to the sum of the values of 424 * its children */ 425 SLN(node_get_mesh(tree, node, &mesh)); 426 ka = sln_mesh_eval(&mesh, nu); 427 if(ka < sum) { 428 ERROR(tree->sln, 429 "ka < ka_{0} + ka_{1} + ... + ka_{N-1}; %e < %e; nu=%-21.20g cm^-1\n", 430 ka, sum, nu); 431 res = RES_BAD_ARG; 432 goto error; 433 } 434 435 /* Complete the probability calculation and normalize the cumulative */ 436 ASSERT(sum != 0); 437 FOR_EACH(i, 0, n) proba[i] /= sum; 438 FOR_EACH(i, 0, n-1) cumul[i] /= sum; 439 cumul[n-1] = 1; /* Handle numerical uncertainty */ 440 441 exit: 442 return res; 443 error: 444 goto exit; 445 } 446 447 static void 448 release_tree(ref_T* ref) 449 { 450 struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref); 451 struct sln_device* sln = NULL; 452 ASSERT(ref); 453 sln = tree->sln; 454 darray_node_release(&tree->nodes); 455 darray_vertex_release(&tree->vertices); 456 if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines)); 457 if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata)); 458 MEM_RM(sln->allocator, tree); 459 SLN(device_ref_put(sln)); 460 } 461 462 /******************************************************************************* 463 * Local function 464 ******************************************************************************/ 465 unsigned 466 node_child_count(const struct sln_node* node, const unsigned tree_arity) 467 { 468 size_t nlines = 0; /* #lines in the node */ 469 size_t nlines_per_child = 0; /* Max #lines in a child */ 470 size_t nchildren = 0; 471 472 /* Pre-conditions */ 473 ASSERT(node && tree_arity >= 2); 474 475 /* Retrieve the node data and compute the #lines it partitions */ 476 nlines = node->range[1] - node->range[0] + 1; 477 ASSERT(nlines); 478 479 /* Based on the arity of the tree, calculate how the lines of the node are 480 * distributed among its children. For low lines count, i.e. when the minimum 481 * number of lines par child is less than the tree arity, the policy below 482 * prioritizes an equal distribution of lines among the children over 483 * maintaining the tree's arity. Thus, if a smaller number of children 484 * results in a more equitable distribution, this option is preferred over 485 * ensuring a number of children equal to the tree's arity. In other words, 486 * the tree's balance is prioritized. */ 487 nlines_per_child = (nlines + tree_arity-1/*ceil*/)/tree_arity; 488 489 /* From the previous line repartition, compute the number of children */ 490 nchildren = (nlines + nlines_per_child-1/*ceil*/)/nlines_per_child; 491 ASSERT(nchildren >= 2); 492 493 ASSERT(nchildren <= UINT_MAX); 494 return (unsigned)nchildren; 495 } 496 497 /******************************************************************************* 498 * Exported symbols 499 ******************************************************************************/ 500 res_T 501 sln_tree_create 502 (struct sln_device* device, 503 const struct sln_tree_create_args* args, 504 struct sln_tree** out_tree) 505 { 506 struct sln_tree* tree = NULL; 507 unsigned nthreads_max = 0; 508 res_T res = RES_OK; 509 510 if(!device || !out_tree) { res = RES_BAD_ARG; goto error; } 511 res = check_sln_tree_create_args(device, FUNC_NAME, args); 512 if(res != RES_OK) goto error; 513 514 res = create_tree(device, FUNC_NAME, &tree); 515 if(res != RES_OK) goto error; 516 SHTR(line_list_ref_get(args->lines)); 517 SHTR(isotope_metadata_ref_get(args->metadata)); 518 tree->args = *args; 519 520 /* Set the #threads to match the maximum number of available threads */ 521 nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs()); 522 tree->args.nthreads_hint = MMIN(tree->args.nthreads_hint, nthreads_max); 523 524 res = tree_build(tree); 525 if(res != RES_OK) goto error; 526 527 exit: 528 if(out_tree) *out_tree = tree; 529 return res; 530 error: 531 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 532 goto exit; 533 } 534 535 res_T 536 sln_tree_read 537 (struct sln_device* sln, 538 const struct sln_tree_read_args* args, 539 struct sln_tree** out_tree) 540 { 541 hash256_T hash_mdata1; 542 hash256_T hash_mdata2; 543 hash256_T hash_lines1; 544 hash256_T hash_lines2; 545 546 struct stream stream = STREAM_NULL; 547 struct sln_tree* tree = NULL; 548 size_t n = 0; 549 int version = 0; 550 res_T res = RES_OK; 551 552 if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; } 553 res = check_sln_tree_read_args(sln, FUNC_NAME, args); 554 if(res != RES_OK) goto error; 555 556 res = create_tree(sln, FUNC_NAME, &tree); 557 if(res != RES_OK) goto error; 558 559 res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream); 560 if(res != RES_OK) goto error; 561 562 #define READ(Var, Nb) { \ 563 if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 564 if(feof(stream.fp)) { \ 565 res = RES_BAD_ARG; \ 566 } else if(ferror(stream.fp)) { \ 567 res = RES_IO_ERR; \ 568 } else { \ 569 res = RES_UNKNOWN_ERR; \ 570 } \ 571 ERROR(sln, "%s: error loading the tree structure -- %s.\n", \ 572 stream.name, res_to_cstr(res)); \ 573 goto error; \ 574 } \ 575 } (void)0 576 READ(&version, 1); 577 if(version != SLN_TREE_VERSION) { 578 ERROR(sln, 579 "%s: unexpected tree version %d. Expecting a tree in version %d.\n", 580 stream.name, version, SLN_TREE_VERSION); 581 res = RES_BAD_ARG; 582 goto error; 583 } 584 585 res = shtr_isotope_metadata_hash(args->metadata, hash_mdata1); 586 if(res != RES_OK) goto error; 587 588 READ(hash_mdata2, sizeof(hash256_T)); 589 if(!hash256_eq(hash_mdata1, hash_mdata2)) { 590 ERROR(sln, 591 "%s: the input isotopic metadata are not those used " 592 "during tree construction.\n", stream.name); 593 res = RES_BAD_ARG; 594 goto error; 595 } 596 597 SHTR(isotope_metadata_ref_get(args->metadata)); 598 tree->args.metadata = args->metadata; 599 600 READ(hash_lines1, sizeof(hash256_T)); 601 if(!args->disable_line_hash_check) { 602 res = shtr_line_list_hash(args->lines, hash_lines2); 603 if(res != RES_OK) goto error; 604 605 if(!hash256_eq(hash_lines1, hash_lines2)) { 606 ERROR(sln, 607 "%s: the input list of lines is not the one used to build the tree.\n", 608 stream.name); 609 res = RES_BAD_ARG; 610 goto error; 611 } 612 } 613 614 SHTR(line_list_ref_get(args->lines)); 615 tree->args.lines = args->lines; 616 617 READ(&n, 1); 618 if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error; 619 READ(darray_node_data_get(&tree->nodes), n); 620 621 READ(&n, 1); 622 if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error; 623 READ(darray_vertex_data_get(&tree->vertices), n); 624 625 READ(&tree->args.line_profile, 1); 626 READ(&tree->args.molecules, 1); 627 READ(&tree->args.pressure, 1); 628 READ(&tree->args.temperature, 1); 629 READ(&tree->args.nvertices_hint, 1); 630 READ(&tree->args.mesh_decimation_err, 1); 631 READ(&tree->args.mesh_type, 1); 632 READ(&tree->args.arity, 1); 633 READ(&tree->args.leaf_nlines, 1); 634 #undef READ 635 636 exit: 637 stream_release(&stream); 638 if(out_tree) *out_tree = tree; 639 return res; 640 error: 641 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 642 goto exit; 643 } 644 645 res_T 646 sln_tree_ref_get(struct sln_tree* tree) 647 { 648 if(!tree) return RES_BAD_ARG; 649 ref_get(&tree->ref); 650 return RES_OK; 651 } 652 653 res_T 654 sln_tree_ref_put(struct sln_tree* tree) 655 { 656 if(!tree) return RES_BAD_ARG; 657 ref_put(&tree->ref, release_tree); 658 return RES_OK; 659 } 660 661 res_T 662 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc) 663 { 664 const struct sln_node* node = NULL; 665 unsigned depth = 0; 666 667 if(!tree || !desc) return RES_BAD_ARG; 668 669 desc->mesh_decimation_err = tree->args.mesh_decimation_err; 670 desc->mesh_type = tree->args.mesh_type; 671 desc->line_profile = tree->args.line_profile; 672 desc->nnodes = darray_node_size_get(&tree->nodes); 673 desc->nvertices = darray_vertex_size_get(&tree->vertices); 674 desc->pressure = tree->args.pressure; /* [atm] */ 675 desc->temperature = tree->args.temperature; /* [K] */ 676 desc->arity = tree->args.arity; 677 desc->leaf_nlines = tree->args.leaf_nlines; 678 679 SHTR(line_list_get_size(tree->args.lines, &desc->nlines)); 680 681 node = sln_tree_get_root(tree); 682 while(!sln_node_is_leaf(node)) { 683 node = sln_node_get_child(tree, node, 0); 684 ++depth; 685 } 686 desc->depth = depth; 687 688 return RES_OK; 689 } 690 691 const struct sln_node* 692 sln_tree_get_root(const struct sln_tree* tree) 693 { 694 ASSERT(tree); 695 if(darray_node_size_get(&tree->nodes)) { 696 return darray_node_cdata_get(&tree->nodes); 697 } else { 698 return NULL; 699 } 700 } 701 702 res_T 703 sln_tree_get_line 704 (const struct sln_tree* tree, 705 const size_t iline, 706 struct sln_line* line) 707 { 708 size_t nlines = 0; 709 res_T res = RES_OK; 710 711 if(!tree || !line) { res = RES_BAD_ARG; goto error; } 712 713 SHTR(line_list_get_size(tree->args.lines, &nlines)); 714 if(iline >= nlines) { res = RES_BAD_ARG; goto error; } 715 716 res = line_setup(tree, iline, line); 717 if(res != RES_OK) { 718 ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n", 719 FUNC_NAME, iline, res_to_cstr(res)); 720 goto error; 721 } 722 723 exit: 724 return res; 725 error: 726 goto exit; 727 } 728 729 int 730 sln_node_is_leaf(const struct sln_node* node) 731 { 732 ASSERT(node); 733 return node->offset == 0; 734 } 735 736 unsigned 737 sln_node_get_child_count 738 (const struct sln_tree* tree, 739 const struct sln_node* node) 740 { 741 ASSERT(tree && node); 742 743 if(sln_node_is_leaf(node)) { 744 return 0; /* No child */ 745 } else { 746 return node_child_count(node, tree->args.arity); 747 } 748 } 749 750 const struct sln_node* 751 sln_node_get_child 752 (const struct sln_tree* tree, 753 const struct sln_node* node, 754 const unsigned ichild) 755 { 756 ASSERT(node && ichild < sln_node_get_child_count(tree, node)); 757 ASSERT(!sln_node_is_leaf(node)); 758 (void)tree; 759 return node + node->offset + ichild; 760 } 761 762 res_T 763 sln_node_get_mesh 764 (const struct sln_tree* tree, 765 const struct sln_node* node, 766 struct sln_mesh* mesh) 767 { 768 if(!tree || !node || !mesh) return RES_BAD_ARG; 769 mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex; 770 mesh->nvertices = node->nvertices; 771 return RES_OK; 772 } 773 774 double 775 sln_node_eval 776 (const struct sln_tree* tree, 777 const struct sln_node* node, 778 const double nu) 779 { 780 double ka = 0; 781 size_t iline; 782 ASSERT(tree && node); 783 784 FOR_EACH(iline, node->range[0], node->range[1]+1) { 785 struct sln_line line = SLN_LINE_NULL; 786 res_T res = RES_OK; 787 788 res = line_setup(tree, iline, &line); 789 if(res != RES_OK) { 790 WARN(tree->sln, "%s: could not setup the line %lu-- %s\n", 791 FUNC_NAME, iline, res_to_cstr(res)); 792 continue; 793 } 794 795 ka += sln_line_eval(tree, &line, nu); 796 } 797 return ka; 798 } 799 800 res_T 801 sln_node_get_desc 802 (const struct sln_tree* tree, 803 const struct sln_node* node, 804 struct sln_node_desc* desc) 805 { 806 if(!tree || !node || !desc) return RES_BAD_ARG; 807 desc->ilines[0] = node->range[0]; 808 desc->ilines[1] = node->range[1]; 809 desc->nvertices = node->nvertices; 810 desc->nchildren = sln_node_get_child_count(tree, node); 811 return RES_OK; 812 } 813 814 const struct sln_node* 815 sln_node_sample_leaf 816 (const struct sln_tree* tree, 817 const struct sln_node* root, 818 const double nu, /*[cm^-1]*/ 819 struct ssp_rng* rng, 820 double* out_leaf_proba) /* May be NULL */ 821 { 822 /* Temporary buffers used to store the cumulative of child nodes and their 823 * probability of being sampled based on their importance */ 824 double cumul[SLN_TREE_ARITY_MAX]; 825 double proba[SLN_TREE_ARITY_MAX]; 826 827 const struct sln_node* node = NULL; 828 double leaf_proba = 1; 829 int depth = 0; 830 res_T res = RES_OK; 831 832 if(!tree || !root || !rng) { res = RES_BAD_ARG; goto error; } 833 834 for(depth=0, node=root; !sln_node_is_leaf(node); ++depth) { 835 double r = 0; /* Random number */ 836 unsigned ichild = 0; 837 838 res = build_node_cumulative(tree, node, nu, proba, cumul); 839 if(res != RES_OK) goto error; 840 841 /* Sample a child node based on its importance. Use a simple linear search, 842 * since the tree's arity is small enough that a binary search is not 843 * necessary. FIXME if performance measurements show that this linear search 844 * incurs a significant cost */ 845 r = ssp_rng_canonical(rng); 846 FOR_EACH(ichild, 0, SLN_TREE_ARITY_MAX) { 847 if(r < cumul[ichild]) { 848 leaf_proba *= proba[ichild]; 849 node = sln_node_get_child(tree, node, ichild); 850 break; 851 } 852 } 853 ASSERT(ichild < SLN_TREE_ARITY_MAX); /* A node should have been sampled */ 854 } 855 856 exit: 857 if(out_leaf_proba) *out_leaf_proba = leaf_proba; 858 return node; 859 error: 860 node = NULL; 861 leaf_proba = NaN; 862 goto exit; 863 } 864 865 double 866 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber) 867 { 868 const struct sln_vertex* vtx0 = NULL; 869 const struct sln_vertex* vtx1 = NULL; 870 const float nu = (float)wavenumber; 871 size_t n; /* #vertices */ 872 double u; /* Linear interpolation parameter */ 873 ASSERT(mesh && mesh->nvertices); 874 875 n = mesh->nvertices; 876 877 /* Handle special cases */ 878 if(n == 1) return mesh->vertices[0].ka; 879 if(nu < mesh->vertices[0].wavenumber 880 || nu > mesh->vertices[n-1].wavenumber) { 881 return 0; 882 } 883 if(nu == mesh->vertices[0].wavenumber) return mesh->vertices[0].ka; 884 if(nu == mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka; 885 886 /* Dichotomic search of the mesh vertex whose wavenumber is greater than or 887 * equal to the submitted wavenumber 'nu' */ 888 vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx); 889 vtx0 = vtx1 - 1; 890 ASSERT(vtx1); /* A vertex is necessary found ...*/ 891 ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */ 892 ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber); 893 894 /* Compute the linear interpolation parameter */ 895 u = (wavenumber - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber); 896 u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */ 897 898 if(u == 0) return vtx0->ka; 899 if(u == 1) return vtx1->ka; 900 return u*(vtx1->ka - vtx0->ka) + vtx0->ka; 901 } 902 903 res_T 904 sln_tree_write 905 (const struct sln_tree* tree, 906 const struct sln_tree_write_args* args) 907 { 908 struct stream stream = STREAM_NULL; 909 size_t nnodes, nverts; 910 hash256_T hash_mdata; 911 hash256_T hash_lines; 912 res_T res = RES_OK; 913 914 if(!tree) { res = RES_BAD_ARG; goto error; } 915 res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args); 916 if(res != RES_OK) goto error; 917 918 res = shtr_isotope_metadata_hash(tree->args.metadata, hash_mdata); 919 if(res != RES_OK) goto error; 920 res = shtr_line_list_hash(tree->args.lines, hash_lines); 921 if(res != RES_OK) goto error; 922 923 res = stream_init 924 (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream); 925 if(res != RES_OK) goto error; 926 927 #define WRITE(Var, Nb) { \ 928 if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 929 ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n", \ 930 FUNC_NAME, stream.name, strerror(errno)); \ 931 res = RES_IO_ERR; \ 932 goto error; \ 933 } \ 934 } (void)0 935 WRITE(&SLN_TREE_VERSION, 1); 936 WRITE(hash_mdata, sizeof(hash256_T)); 937 WRITE(hash_lines, sizeof(hash256_T)); 938 939 nnodes = darray_node_size_get(&tree->nodes); 940 WRITE(&nnodes, 1); 941 WRITE(darray_node_cdata_get(&tree->nodes), nnodes); 942 943 nverts = darray_vertex_size_get(&tree->vertices); 944 WRITE(&nverts, 1); 945 WRITE(darray_vertex_cdata_get(&tree->vertices), nverts); 946 947 WRITE(&tree->args.line_profile, 1); 948 WRITE(&tree->args.molecules, 1); 949 WRITE(&tree->args.pressure, 1); 950 WRITE(&tree->args.temperature, 1); 951 WRITE(&tree->args.nvertices_hint, 1); 952 WRITE(&tree->args.mesh_decimation_err, 1); 953 WRITE(&tree->args.mesh_type, 1); 954 WRITE(&tree->args.arity, 1); 955 WRITE(&tree->args.leaf_nlines, 1); 956 #undef WRITE 957 958 exit: 959 stream_release(&stream); 960 return res; 961 error: 962 goto exit; 963 }