test_sln_tree.c (14758B)
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 "test_sln_lines.h" 20 #include "sln.h" 21 22 #include <star/shtr.h> 23 24 #include <rsys/algorithm.h> 25 #include <rsys/math.h> 26 #include <rsys/mem_allocator.h> 27 28 /******************************************************************************* 29 * Helper function 30 ******************************************************************************/ 31 /* Return the index of the line in the line_list or SIZE_MAX if the line does 32 * not exist */ 33 static INLINE size_t 34 find_line 35 (struct shtr_line_list* line_list, 36 const struct shtr_line* line) 37 { 38 struct shtr_line ln = SHTR_LINE_NULL; 39 size_t lo, hi, mid; 40 size_t iline; 41 size_t nlines; 42 43 CHK(shtr_line_list_get_size(line_list, &nlines) == RES_OK); 44 45 /* Dichotomic search */ 46 lo = 0; hi = nlines -1; 47 while(lo < hi) { 48 mid = (lo+hi)/2; 49 50 CHK(shtr_line_list_at(line_list, mid, &ln) == RES_OK); 51 if(line->wavenumber > ln.wavenumber) { 52 lo = mid + 1; 53 } else { 54 hi = mid; 55 } 56 } 57 iline = lo; 58 59 CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK); 60 if(ln.wavenumber != line->wavenumber) return SIZE_MAX; 61 62 63 /* Find a line with the same wavenumber as the one searched for and whose 64 * other member variables are also equal to those of the line searched for */ 65 while(ln.wavenumber == line->wavenumber 66 && !shtr_line_eq(&ln, line) 67 && iline < nlines) { 68 iline += 1; 69 CHK(shtr_line_list_at(line_list, iline, &ln) == RES_OK); 70 } 71 72 return shtr_line_eq(&ln, line) ? iline : SIZE_MAX; 73 } 74 75 /* This test assumes that all the lines contained into the list are 76 * partitionned in the tree */ 77 static void 78 check_tree 79 (const struct sln_tree* tree, 80 struct shtr_line_list* line_list) 81 { 82 #define STACK_SIZE 64 83 const struct sln_node* stack[STACK_SIZE] = {NULL}; 84 struct sln_tree_desc desc = SLN_TREE_DESC_NULL; 85 size_t istack = 0; 86 const struct sln_node* node = NULL; 87 88 char* found_lines = NULL; 89 size_t found_nlines = 0; 90 size_t line_list_sz; 91 92 CHK(shtr_line_list_get_size(line_list, &line_list_sz) == RES_OK); 93 CHK(sln_tree_get_desc(tree, &desc) == RES_OK); 94 95 CHK(found_lines = mem_calloc(line_list_sz, sizeof(char))); 96 97 node = sln_tree_get_root(tree); 98 while(node) { 99 if(!sln_node_is_leaf(node)) { 100 CHK(sln_node_get_lines_count(node) > desc.max_nlines_per_leaf); 101 102 ASSERT(istack < STACK_SIZE); 103 stack[istack++] = sln_node_get_child(node, 1); /* Push the child 1 */ 104 node = sln_node_get_child(node, 0); /* Handle the child 0 */ 105 106 } else { 107 size_t iline, nlines; 108 109 nlines = sln_node_get_lines_count(node); 110 CHK(nlines <= desc.max_nlines_per_leaf); 111 FOR_EACH(iline, 0, nlines) { 112 struct shtr_line line = SHTR_LINE_NULL; 113 size_t found_iline = 0; 114 CHK(sln_node_get_line(tree, node, iline, &line) == RES_OK); 115 found_iline = find_line(line_list, &line); 116 CHK(found_iline != SIZE_MAX); /* Line is found */ 117 118 if(!found_lines[found_iline]) { 119 found_nlines += 1; 120 found_lines[found_iline] = 1; 121 } 122 } 123 124 node = istack ? stack[--istack] : NULL; /* Pop the next node */ 125 } 126 } 127 128 /* Check that all lines are found */ 129 CHK(found_nlines == line_list_sz); 130 131 mem_rm(found_lines); 132 #undef STACK_SIZE 133 } 134 135 static void 136 dump_line(FILE* stream, const struct sln_mesh* mesh) 137 { 138 size_t i; 139 CHK(mesh); 140 FOR_EACH(i, 0, mesh->nvertices) { 141 fprintf(stream, "%g %g\n", 142 mesh->vertices[i].wavenumber, 143 mesh->vertices[i].ka); 144 } 145 } 146 147 static void 148 test_tree 149 (struct sln_device* sln, 150 const struct sln_tree_create_args* tree_args_in, 151 struct shtr_line_list* line_list) 152 { 153 struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT; 154 struct sln_tree_desc desc = SLN_TREE_DESC_NULL; 155 struct sln_node_desc node_desc = SLN_NODE_DESC_NULL; 156 struct sln_mesh mesh = SLN_MESH_NULL; 157 struct sln_tree* tree = NULL; 158 const struct sln_node* node = NULL; 159 struct shtr_line line = SHTR_LINE_NULL; 160 size_t nlines = 0; 161 162 CHK(sln && tree_args_in && line_list); 163 tree_args = *tree_args_in; 164 165 CHK(sln_tree_create(NULL, &tree_args, &tree) == RES_BAD_ARG); 166 CHK(sln_tree_create(sln, NULL, &tree) == RES_BAD_ARG); 167 CHK(sln_tree_create(sln, &tree_args, NULL) == RES_BAD_ARG); 168 CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK); 169 170 CHK(shtr_line_list_get_size(line_list, &nlines) == RES_OK); 171 172 CHK(sln_tree_get_desc(NULL, &desc) == RES_BAD_ARG); 173 CHK(sln_tree_get_desc(tree, NULL) == RES_BAD_ARG); 174 CHK(sln_tree_get_desc(tree, &desc) == RES_OK); 175 176 CHK(desc.max_nlines_per_leaf >= 1); 177 CHK(desc.mesh_decimation_err == tree_args.mesh_decimation_err); 178 CHK(desc.mesh_type == tree_args.mesh_type); 179 CHK(desc.line_profile == tree_args.line_profile); 180 CHK(desc.nlines == nlines); 181 CHK(desc.nvertices >= 1); 182 CHK(desc.nnodes >= 1); 183 CHK(desc.depth == (unsigned)log2i((int)round_up_pow2(desc.nlines))); 184 185 CHK(node = sln_tree_get_root(tree)); 186 CHK(node != NULL); 187 188 CHK(sln_node_get_desc(NULL, &node_desc) == RES_BAD_ARG); 189 CHK(sln_node_get_desc(node, NULL) == RES_BAD_ARG); 190 CHK(sln_node_get_desc(node, &node_desc) == RES_OK); 191 CHK(node_desc.nlines = desc.nlines); 192 CHK(node_desc.nvertices >= 1 && node_desc.nvertices < desc.nvertices); 193 194 while(!sln_node_is_leaf(node)) { 195 struct sln_node_desc node_desc_next = SLN_NODE_DESC_NULL; 196 node = sln_node_get_child(node, 0); 197 198 CHK(sln_node_get_desc(node, &node_desc_next) == RES_OK); 199 CHK(node_desc_next.nlines >= 1); 200 CHK(node_desc_next.nlines < node_desc.nlines); 201 CHK(node_desc_next.nvertices >= 1); 202 CHK(node_desc_next.nlines < node_desc.nvertices); 203 204 node_desc = node_desc_next; 205 } 206 207 CHK(sln_node_get_lines_count(node) <= desc.max_nlines_per_leaf); 208 CHK(sln_node_get_line(NULL, node, 0, &line) == RES_BAD_ARG); 209 CHK(sln_node_get_line(tree, NULL, 0, &line) == RES_BAD_ARG); 210 CHK(sln_node_get_line(tree, node, desc.max_nlines_per_leaf+1, &line) 211 == RES_BAD_ARG); 212 CHK(sln_node_get_line(tree, node, 0, NULL) == RES_BAD_ARG); 213 CHK(sln_node_get_line(tree, node, 0, &line) == RES_OK); 214 CHK(sln_node_get_line(tree, sln_tree_get_root(tree), 0, &line) == RES_OK); 215 216 CHK(sln_node_get_mesh(NULL, node, &mesh) == RES_BAD_ARG); 217 CHK(sln_node_get_mesh(tree, NULL, &mesh) == RES_BAD_ARG); 218 CHK(sln_node_get_mesh(tree, node, NULL) == RES_BAD_ARG); 219 CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK); 220 221 CHK(node = sln_tree_get_root(tree)); 222 CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK); 223 224 dump_line(stdout, &mesh); 225 check_tree(tree, line_list); 226 227 CHK(sln_tree_ref_get(NULL) == RES_BAD_ARG); 228 CHK(sln_tree_ref_get(tree) == RES_OK); 229 CHK(sln_tree_ref_put(NULL) == RES_BAD_ARG); 230 CHK(sln_tree_ref_put(tree) == RES_OK); 231 CHK(sln_tree_ref_put(tree) == RES_OK); 232 233 tree_args.mesh_decimation_err = -1; 234 CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG); 235 236 tree_args.mesh_decimation_err = 1e-1f; 237 tree_args.mesh_type = SLN_MESH_TYPES_COUNT__; 238 CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG); 239 240 tree_args.mesh_type = SLN_MESH_FIT; 241 tree_args.line_profile = SLN_LINE_PROFILES_COUNT__; 242 CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG); 243 244 tree_args.line_profile = SLN_LINE_PROFILE_VOIGT; 245 CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK); 246 CHK(sln_tree_ref_put(tree) == RES_OK); 247 } 248 249 static void 250 check_mesh_equality(const struct sln_mesh* mesh1, const struct sln_mesh* mesh2) 251 { 252 size_t i; 253 CHK(mesh1 && mesh2); 254 CHK(mesh1->nvertices == mesh2->nvertices); 255 256 FOR_EACH(i, 0, mesh1->nvertices) { 257 CHK(mesh1->vertices[i].wavenumber == mesh2->vertices[i].wavenumber); 258 CHK(mesh1->vertices[i].ka == mesh2->vertices[i].ka); 259 } 260 } 261 262 static void 263 check_node_equality 264 (const struct sln_tree* tree1, 265 const struct sln_tree* tree2, 266 const struct sln_node* node1, 267 const struct sln_node* node2) 268 { 269 struct sln_mesh mesh1 = SLN_MESH_NULL; 270 struct sln_mesh mesh2 = SLN_MESH_NULL; 271 size_t iline, nlines; 272 273 CHK(node1 && node2); 274 CHK(sln_node_is_leaf(node1) == sln_node_is_leaf(node2)); 275 CHK(sln_node_get_lines_count(node1) == sln_node_get_lines_count(node2)); 276 nlines = sln_node_get_lines_count(node1); 277 278 FOR_EACH(iline, 0, nlines) { 279 struct shtr_line line1 = SHTR_LINE_NULL; 280 struct shtr_line line2 = SHTR_LINE_NULL; 281 CHK(sln_node_get_line(tree1, node1, iline, &line1) == RES_OK); 282 CHK(sln_node_get_line(tree2, node2, iline, &line2) == RES_OK); 283 CHK(shtr_line_eq(&line1, &line2)); 284 } 285 286 CHK(sln_node_get_mesh(tree1, node1, &mesh1) == RES_OK); 287 CHK(sln_node_get_mesh(tree2, node2, &mesh2) == RES_OK); 288 check_mesh_equality(&mesh1, &mesh2); 289 } 290 291 static void 292 check_tree_equality 293 (const struct sln_tree* tree1, 294 const struct sln_tree* tree2) 295 { 296 const struct sln_node* stack[128] = {NULL}; 297 int istack = 0; 298 299 struct sln_tree_desc desc1 = SLN_TREE_DESC_NULL; 300 struct sln_tree_desc desc2 = SLN_TREE_DESC_NULL; 301 const struct sln_node* node1 = NULL; 302 const struct sln_node* node2 = NULL; 303 304 CHK(sln_tree_get_desc(tree1, &desc1) == RES_OK); 305 CHK(sln_tree_get_desc(tree2, &desc2) == RES_OK); 306 CHK(desc1.max_nlines_per_leaf == desc2.max_nlines_per_leaf); 307 CHK(desc1.mesh_decimation_err == desc2.mesh_decimation_err); 308 CHK(desc1.line_profile == desc2.line_profile); 309 310 stack[istack++] = NULL; 311 stack[istack++] = NULL; 312 313 node1 = sln_tree_get_root(tree1); 314 node2 = sln_tree_get_root(tree2); 315 316 while(node1 || node2) { 317 CHK((!node1 && !node2) || (node1 && node2)); 318 check_node_equality(tree1, tree2, node1, node2); 319 320 if(sln_node_is_leaf(node1)) { 321 node2 = stack[--istack]; 322 node1 = stack[--istack]; 323 } else { 324 stack[istack++] = sln_node_get_child(node1, 1); 325 stack[istack++] = sln_node_get_child(node2, 1); 326 node1 = sln_node_get_child(node1, 0); 327 node2 = sln_node_get_child(node2, 0); 328 } 329 } 330 } 331 332 static void 333 test_tree_serialization 334 (struct sln_device* sln, 335 const struct sln_tree_create_args* tree_args, 336 struct shtr* shtr) 337 { 338 struct sln_tree_write_args wargs = SLN_TREE_WRITE_ARGS_NULL; 339 struct sln_tree_read_args rargs = SLN_TREE_READ_ARGS_NULL; 340 struct sln_tree* tree1 = NULL; 341 struct sln_tree* tree2 = NULL; 342 343 const char* filename = "tree.sln"; 344 FILE* fp = NULL; 345 346 CHK(sln_tree_create(sln, tree_args, &tree1) == RES_OK); 347 348 CHK(fp = fopen(filename, "w+")); 349 350 wargs.file = fp; 351 CHK(sln_tree_write(NULL, &wargs) == RES_BAD_ARG); 352 CHK(sln_tree_write(tree1, NULL) == RES_BAD_ARG); 353 wargs.file = NULL; 354 CHK(sln_tree_write(tree1, &wargs) == RES_BAD_ARG); 355 wargs.file = fp; 356 CHK(sln_tree_write(tree1, &wargs) == RES_OK); 357 rewind(fp); 358 359 rargs.shtr = shtr; 360 rargs.file = fp; 361 CHK(sln_tree_read(NULL, &rargs, &tree2) == RES_BAD_ARG); 362 CHK(sln_tree_read(sln, NULL, &tree2) == RES_BAD_ARG); 363 rargs.shtr = NULL; 364 CHK(sln_tree_read(sln, &rargs, &tree2) == RES_BAD_ARG); 365 rargs.shtr = shtr; 366 rargs.file = NULL; 367 CHK(sln_tree_read(sln, &rargs, &tree2) == RES_BAD_ARG); 368 rargs.file = fp; 369 CHK(sln_tree_read(sln, &rargs, NULL) == RES_BAD_ARG); 370 CHK(sln_tree_read(sln, &rargs, &tree2) == RES_OK); 371 fclose(fp); 372 373 check_tree_equality(tree1, tree2); 374 CHK(sln_tree_ref_put(tree2) == RES_OK); 375 376 wargs.file = NULL; 377 wargs.filename = filename; 378 CHK(sln_tree_write(tree1, &wargs) == RES_OK); 379 380 rargs.file = NULL; 381 rargs.filename = "nop"; 382 CHK(sln_tree_read(sln, &rargs, &tree2) == RES_IO_ERR); 383 rargs.filename = filename; 384 CHK(sln_tree_read(sln, &rargs, &tree2) == RES_OK); 385 386 check_tree_equality(tree1, tree2); 387 388 CHK(sln_tree_ref_put(tree2) == RES_OK); 389 CHK(sln_tree_ref_put(tree1) == RES_OK); 390 } 391 392 /******************************************************************************* 393 * Test function 394 ******************************************************************************/ 395 int 396 main(int argc, char** argv) 397 { 398 struct sln_device_create_args dev_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 399 struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT; 400 struct sln_device* sln = NULL; 401 402 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 403 struct shtr* shtr = NULL; 404 struct shtr_line_list_load_args line_load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL; 405 struct shtr_line_list* line_list = NULL; 406 struct shtr_isotope_metadata* metadata = NULL; 407 408 FILE* fp_lines = NULL; 409 FILE* fp_mdata = NULL; 410 (void)argc, (void)argv; 411 412 /* Generate the file of the isotope metadata */ 413 CHK(fp_mdata = tmpfile()); 414 fprintf(fp_mdata, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); 415 write_shtr_molecule(fp_mdata, &g_H2O); 416 write_shtr_molecule(fp_mdata, &g_CO2); 417 write_shtr_molecule(fp_mdata, &g_O3); 418 rewind(fp_mdata); 419 420 /* Generate the file of lines */ 421 CHK(fp_lines = tmpfile()); 422 write_shtr_lines(fp_lines, g_lines, g_nlines); 423 rewind(fp_lines); 424 425 /* Load the isotope metadata and the lines */ 426 shtr_args.verbose = 1; 427 CHK(shtr_create(&shtr_args, &shtr) == RES_OK); 428 CHK(shtr_isotope_metadata_load_stream(shtr, fp_mdata, NULL, &metadata) == RES_OK); 429 430 line_load_args.filename = "stream"; 431 line_load_args.file = fp_lines; 432 CHK(shtr_line_list_load(shtr, &line_load_args, &line_list) == RES_OK); 433 434 CHK(fclose(fp_lines) == 0); 435 CHK(fclose(fp_mdata) == 0); 436 437 dev_args.verbose = 1; 438 CHK(sln_device_create(&dev_args, &sln) == RES_OK); 439 440 /* Create the mixture */ 441 tree_args.metadata = metadata; 442 tree_args.lines = line_list; 443 tree_args.molecules[SHTR_H2O].concentration = 1.0/3.0; 444 tree_args.molecules[SHTR_H2O].cutoff = 25; 445 tree_args.molecules[SHTR_CO2].concentration = 1.0/3.0; 446 tree_args.molecules[SHTR_CO2].cutoff = 50; 447 tree_args.molecules[SHTR_O3 ].concentration = 1.0/3.0; 448 tree_args.molecules[SHTR_O3 ].cutoff = 25; 449 tree_args.pressure = 1; 450 tree_args.temperature = 296; 451 452 test_tree(sln, &tree_args, line_list); 453 test_tree_serialization(sln, &tree_args, shtr); 454 455 CHK(sln_device_ref_put(sln) == RES_OK); 456 CHK(shtr_ref_put(shtr) == RES_OK); 457 CHK(shtr_line_list_ref_put(line_list) == RES_OK); 458 CHK(shtr_isotope_metadata_ref_put(metadata) == RES_OK); 459 CHK(mem_allocated_size() == 0); 460 return 0; 461 }