sln_tree.c (16085B)
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 26 #include <rsys/algorithm.h> 27 #include <rsys/cstr.h> 28 #include <rsys/math.h> 29 30 /******************************************************************************* 31 * Helper functions 32 ******************************************************************************/ 33 /* Check that the sum of the molecular concentrations is equal to 1 */ 34 static res_T 35 check_molecule_concentration 36 (struct sln_device* sln, 37 const char* caller, 38 const struct sln_tree_create_args* args) 39 { 40 int i = 0; 41 double sum = 0; 42 ASSERT(sln && caller && args); 43 44 FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) { 45 if(i == SHTR_MOLECULE_ID_NULL) continue; 46 sum += args->molecules[i].concentration; 47 } 48 49 if(!eq_eps(sum, 1, 1e-6)) { 50 ERROR(sln, 51 "%s: the sum of the concentrations of the molecules does not equal 1: %g.\n", 52 caller, sum); 53 return RES_BAD_ARG; 54 } 55 56 return RES_OK; 57 } 58 59 /* Verify that the isotope abundance are valids */ 60 static res_T 61 check_molecule_isotope_abundances 62 (struct sln_device* sln, 63 const char* caller, 64 const struct sln_molecule* molecule) 65 { 66 int i = 0; 67 double sum = 0; 68 ASSERT(sln && caller && molecule); 69 70 /* The isotopic abundances are the default ones. Nothing to do */ 71 if(!molecule->non_default_isotope_abundances) return RES_OK; 72 73 /* The isotopic abundances are not the default ones. 74 * Verify that they are valid ... */ 75 FOR_EACH(i, 0, SLN_MAX_ISOTOPES_COUNT) { 76 if(molecule->isotope_abundances[i] < 0) { 77 const int isotope_id = i + 1; /* isotope id in [1, 10] */ 78 ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n", 79 caller, isotope_id, shtr_molecule_cstr(i), 80 molecule->isotope_abundances[i]); 81 return RES_BAD_ARG; 82 } 83 84 sum += molecule->isotope_abundances[i]; 85 } 86 87 /* ... and that their sum equals 1 */ 88 if(eq_eps(sum, 1, 1e-6)) { 89 ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n", 90 caller, shtr_molecule_cstr(i), sum); 91 return RES_BAD_ARG; 92 } 93 94 return RES_OK; 95 } 96 97 static res_T 98 check_molecules 99 (struct sln_device* sln, 100 const char* caller, 101 const struct sln_tree_create_args* args) 102 { 103 char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0}; 104 105 size_t iline = 0; 106 size_t nlines = 0; 107 res_T res = RES_OK; 108 ASSERT(args->lines); 109 110 res = check_molecule_concentration(sln, caller, args); 111 if(res != RES_OK) return res; 112 113 /* Interate over the lines to define which molecules has to be checked, i.e., 114 * the ones used in the mixture */ 115 SHTR(line_list_get_size(args->lines, &nlines)); 116 FOR_EACH(iline, 0, nlines) { 117 struct shtr_line line = SHTR_LINE_NULL; 118 const struct sln_molecule* molecule = NULL; 119 120 SHTR(line_list_at(args->lines, iline, &line)); 121 122 /* This molecule was already checked */ 123 if(molecule_ok[line.molecule_id]) continue; 124 125 molecule = args->molecules + line.molecule_id; 126 127 if(molecule->concentration == 0) { 128 /* A molecular concentration of zero is allowed, but may be a user error, 129 * as 0 is the default concentration in the tree creation arguments. 130 * Therefore, warn the user about this value so that they can determine 131 * whether or not it is an error on their part. */ 132 WARN(sln, "%s: the concentration of %s is zero\n.\n", 133 caller, shtr_molecule_cstr(line.molecule_id)); 134 135 } else if(molecule->concentration < 0) { 136 /* Concentration cannot be negative... */ 137 ERROR(sln, "%s: invalid %s concentration: %g.\n", 138 FUNC_NAME, shtr_molecule_cstr(line.molecule_id), 139 molecule->concentration); 140 return RES_BAD_ARG; 141 } 142 143 if(molecule->cutoff <= 0) { 144 /* ... cutoff either */ 145 ERROR(sln, "%s: invalid %s cutoff: %g.\n", 146 caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff); 147 return RES_BAD_ARG; 148 } 149 150 res = check_molecule_isotope_abundances(sln, caller, molecule); 151 if(res != RES_OK) return res; 152 153 molecule_ok[line.molecule_id] = 1; 154 } 155 156 return RES_OK; 157 } 158 159 static res_T 160 check_sln_tree_create_args 161 (struct sln_device* sln, 162 const char* caller, 163 const struct sln_tree_create_args* args) 164 { 165 res_T res = RES_OK; 166 ASSERT(sln && caller); 167 168 if(!args) return RES_BAD_ARG; 169 170 if(!args->metadata) { 171 ERROR(sln, "%s: the metadata is missing.\n", caller); 172 return RES_BAD_ARG; 173 } 174 175 if(!args->lines) { 176 ERROR(sln, "%s: the list of lines is missing.\n", caller); 177 return RES_BAD_ARG; 178 } 179 180 if(args->nvertices_hint == 0) { 181 ERROR(sln, 182 "%s: invalid hint on the number of vertices around the line center %lu.\n", 183 caller, (unsigned long)args->nvertices_hint); 184 return RES_BAD_ARG; 185 } 186 187 if(args->mesh_decimation_err < 0) { 188 ERROR(sln, "%s: invalid decimation error %g.\n", 189 caller, args->mesh_decimation_err); 190 return RES_BAD_ARG; 191 } 192 193 if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) { 194 ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type); 195 return RES_BAD_ARG; 196 } 197 198 if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) { 199 ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile); 200 return RES_BAD_ARG; 201 } 202 203 res = check_molecules(sln, caller, args); 204 if(res != RES_OK) return res; 205 206 return RES_OK; 207 } 208 209 static res_T 210 create_tree 211 (struct sln_device* sln, 212 const char* caller, 213 struct sln_tree** out_tree) 214 { 215 struct sln_tree* tree = NULL; 216 res_T res = RES_OK; 217 ASSERT(sln && caller && out_tree); 218 219 tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree)); 220 if(!tree) { 221 ERROR(sln, "%s: could not allocate the tree data structure.\n", 222 caller); 223 res = RES_MEM_ERR; 224 goto error; 225 } 226 ref_init(&tree->ref); 227 SLN(device_ref_get(sln)); 228 tree->sln = sln; 229 darray_node_init(sln->allocator, &tree->nodes); 230 darray_vertex_init(sln->allocator, &tree->vertices); 231 232 exit: 233 *out_tree = tree; 234 return res; 235 error: 236 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 237 goto exit; 238 } 239 240 static INLINE int 241 cmp_nu_vtx(const void* key, const void* item) 242 { 243 const float nu = *((const float*)key); 244 const struct sln_vertex* vtx = item; 245 if(nu < vtx->wavenumber) return -1; 246 if(nu > vtx->wavenumber) return +1; 247 return 0; 248 } 249 250 static void 251 release_tree(ref_T* ref) 252 { 253 struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref); 254 struct sln_device* sln = NULL; 255 ASSERT(ref); 256 sln = tree->sln; 257 darray_node_release(&tree->nodes); 258 darray_vertex_release(&tree->vertices); 259 if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines)); 260 if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata)); 261 MEM_RM(sln->allocator, tree); 262 SLN(device_ref_put(sln)); 263 } 264 265 /******************************************************************************* 266 * Exported symbols 267 ******************************************************************************/ 268 res_T 269 sln_tree_create 270 (struct sln_device* device, 271 const struct sln_tree_create_args* args, 272 struct sln_tree** out_tree) 273 { 274 struct sln_tree* tree = NULL; 275 res_T res = RES_OK; 276 277 if(!device || !out_tree) { res = RES_BAD_ARG; goto error; } 278 res = check_sln_tree_create_args(device, FUNC_NAME, args); 279 if(res != RES_OK) goto error; 280 281 res = create_tree(device, FUNC_NAME, &tree); 282 if(res != RES_OK) goto error; 283 SHTR(line_list_ref_get(args->lines)); 284 SHTR(isotope_metadata_ref_get(args->metadata)); 285 tree->args = *args; 286 287 res = tree_build(tree); 288 if(res != RES_OK) goto error; 289 290 exit: 291 if(out_tree) *out_tree = tree; 292 return res; 293 error: 294 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 295 goto exit; 296 } 297 298 res_T 299 sln_tree_create_from_stream 300 (struct sln_device* sln, 301 struct shtr* shtr, 302 FILE* stream, 303 struct sln_tree** out_tree) 304 { 305 struct sln_tree* tree = NULL; 306 size_t n = 0; 307 int version = 0; 308 res_T res = RES_OK; 309 310 if(!sln || !shtr || !stream || !out_tree) { 311 res = RES_BAD_ARG; 312 goto error; 313 } 314 315 res = create_tree(sln, FUNC_NAME, &tree); 316 if(res != RES_OK) goto error; 317 318 #define READ(Var, Nb) { \ 319 if(fread((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \ 320 if(feof(stream)) { \ 321 res = RES_BAD_ARG; \ 322 } else if(ferror(stream)) { \ 323 res = RES_IO_ERR; \ 324 } else { \ 325 res = RES_UNKNOWN_ERR; \ 326 } \ 327 goto error; \ 328 } \ 329 } (void)0 330 READ(&version, 1); 331 if(version != SLN_TREE_VERSION) { 332 ERROR(sln, 333 "%s: unexpected tree version %d. Expecting a tree in version %d.\n", 334 FUNC_NAME, version, SLN_TREE_VERSION); 335 res = RES_BAD_ARG; 336 goto error; 337 } 338 339 READ(&n, 1); 340 if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error; 341 READ(darray_node_data_get(&tree->nodes), n); 342 343 READ(&n, 1); 344 if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error; 345 READ(darray_vertex_data_get(&tree->vertices), n); 346 347 READ(&tree->args.line_profile, 1); 348 READ(&tree->args.molecules, 1); 349 READ(&tree->args.pressure, 1); 350 READ(&tree->args.temperature, 1); 351 READ(&tree->args.nvertices_hint, 1); 352 READ(&tree->args.mesh_decimation_err, 1); 353 READ(&tree->args.mesh_type, 1); 354 #undef READ 355 356 res = shtr_isotope_metadata_create_from_stream(shtr, stream, &tree->args.metadata); 357 if(res != RES_OK) goto error; 358 359 res = shtr_line_list_create_from_stream(shtr, stream, &tree->args.lines); 360 if(res != RES_OK) goto error; 361 362 exit: 363 if(out_tree) *out_tree = tree; 364 return res; 365 error: 366 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 367 if(sln) { 368 ERROR(sln, "%s: error loading the tree structure -- %s.\n", 369 FUNC_NAME, res_to_cstr(res)); 370 } 371 goto exit; 372 } 373 374 res_T 375 sln_tree_ref_get(struct sln_tree* tree) 376 { 377 if(!tree) return RES_BAD_ARG; 378 ref_get(&tree->ref); 379 return RES_OK; 380 } 381 382 res_T 383 sln_tree_ref_put(struct sln_tree* tree) 384 { 385 if(!tree) return RES_BAD_ARG; 386 ref_put(&tree->ref, release_tree); 387 return RES_OK; 388 } 389 390 res_T 391 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc) 392 { 393 if(!tree || !desc) return RES_BAD_ARG; 394 desc->max_nlines_per_leaf = 1; 395 desc->mesh_decimation_err = tree->args.mesh_decimation_err; 396 desc->mesh_type = tree->args.mesh_type; 397 desc->line_profile = tree->args.line_profile; 398 return RES_OK; 399 } 400 401 const struct sln_node* 402 sln_tree_get_root(const struct sln_tree* tree) 403 { 404 ASSERT(tree); 405 if(darray_node_size_get(&tree->nodes)) { 406 return darray_node_cdata_get(&tree->nodes); 407 } else { 408 return NULL; 409 } 410 } 411 412 int 413 sln_node_is_leaf(const struct sln_node* node) 414 { 415 ASSERT(node); 416 return node->offset == 0; 417 } 418 419 const struct sln_node* 420 sln_node_get_child(const struct sln_node* node, const unsigned ichild) 421 { 422 ASSERT(node && ichild <= 1); 423 return node + node->offset + ichild; 424 } 425 426 size_t 427 sln_node_get_lines_count(const struct sln_node* node) 428 { 429 ASSERT(node); 430 return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/; 431 } 432 433 res_T 434 sln_node_get_line 435 (const struct sln_tree* tree, 436 const struct sln_node* node, 437 const size_t iline, 438 struct shtr_line* line) 439 { 440 if(!tree || !node || iline > sln_node_get_lines_count(node)) 441 return RES_BAD_ARG; 442 443 return shtr_line_list_at 444 (tree->args.lines, node->range[0] + iline, line); 445 } 446 447 res_T 448 sln_node_get_mesh 449 (const struct sln_tree* tree, 450 const struct sln_node* node, 451 struct sln_mesh* mesh) 452 { 453 if(!tree || !node || !mesh) return RES_BAD_ARG; 454 mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex; 455 mesh->nvertices = node->nvertices; 456 return RES_OK; 457 } 458 459 double 460 sln_node_eval 461 (const struct sln_tree* tree, 462 const struct sln_node* node, 463 const double nu) 464 { 465 double ka = 0; 466 size_t iline; 467 ASSERT(tree && node); 468 469 FOR_EACH(iline, node->range[0], node->range[1]+1) { 470 struct line line = LINE_NULL; 471 res_T res = RES_OK; 472 473 res = line_setup(tree, iline, &line); 474 if(res != RES_OK) { 475 WARN(tree->sln, "%s: could not setup the line %lu-- %s\n", 476 FUNC_NAME, iline, res_to_cstr(res)); 477 } 478 479 ka += line_value(tree, &line, nu); 480 } 481 return ka; 482 } 483 484 double 485 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber) 486 { 487 const struct sln_vertex* vtx0 = NULL; 488 const struct sln_vertex* vtx1 = NULL; 489 const float nu = (float)wavenumber; 490 size_t n; /* #vertices */ 491 float u; /* Linear interpolation paraeter */ 492 ASSERT(mesh && mesh->nvertices); 493 494 n = mesh->nvertices; 495 496 /* Handle special cases */ 497 if(n == 1) return mesh->vertices[0].ka; 498 if(nu <= mesh->vertices[0].wavenumber) return mesh->vertices[0].ka; 499 if(nu >= mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka; 500 501 /* Dichotomic search of the mesh vertex whose wavenumber is greater than or 502 * equal to the submitted wavenumber 'nu' */ 503 vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx); 504 vtx0 = vtx1 - 1; 505 ASSERT(vtx1); /* A vertex is necessary found ...*/ 506 ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */ 507 ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber); 508 509 /* Compute the linear interpolation parameter */ 510 u = (nu - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber); 511 u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */ 512 513 if(u == 0) return vtx0->ka; 514 if(u == 1) return vtx1->ka; 515 return u*(vtx1->ka - vtx0->ka) + vtx0->ka; 516 } 517 518 res_T 519 sln_tree_write(const struct sln_tree* tree, FILE* stream) 520 { 521 size_t n; 522 res_T res = RES_OK; 523 524 if(!tree || !stream) { res = RES_BAD_ARG; goto error; } 525 526 #define WRITE(Var, Nb) { \ 527 if(fwrite((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \ 528 res = RES_IO_ERR; \ 529 goto error; \ 530 } \ 531 } (void)0 532 WRITE(&SLN_TREE_VERSION, 1); 533 534 n = darray_node_size_get(&tree->nodes); 535 WRITE(&n, 1); 536 WRITE(darray_node_cdata_get(&tree->nodes), n); 537 538 n = darray_vertex_size_get(&tree->vertices); 539 WRITE(&n, 1); 540 WRITE(darray_vertex_cdata_get(&tree->vertices), n); 541 542 WRITE(&tree->args.line_profile, 1); 543 WRITE(&tree->args.molecules, 1); 544 WRITE(&tree->args.pressure, 1); 545 WRITE(&tree->args.temperature, 1); 546 WRITE(&tree->args.nvertices_hint, 1); 547 WRITE(&tree->args.mesh_decimation_err, 1); 548 WRITE(&tree->args.mesh_type, 1); 549 #undef WRITE 550 551 res = shtr_isotope_metadata_write(tree->args.metadata, stream); 552 if(res != RES_OK) goto error; 553 554 res = shtr_line_list_write(tree->args.lines, stream); 555 if(res != RES_OK) goto error; 556 557 exit: 558 return res; 559 error: 560 if(tree) { 561 ERROR(tree->sln, "%s: error writing the tree -- %s\n", 562 FUNC_NAME, res_to_cstr(res)); 563 } 564 goto exit; 565 }