sln_tree_build.c (16305B)
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_polyline.h" 22 #include "sln_tree_c.h" 23 24 #include <star/shtr.h> 25 26 #include <rsys/cstr.h> 27 28 /******************************************************************************* 29 * Helper functions 30 ******************************************************************************/ 31 static FINLINE uint32_t 32 ui64_to_ui32(const uint64_t ui64) 33 { 34 if(ui64 > UINT32_MAX) 35 VFATAL("%s: overflow %lu.\n", ARG2(FUNC_NAME, ((unsigned long)ui64))); 36 return (uint32_t)ui64; 37 } 38 39 static INLINE res_T 40 build_leaf_polyline(struct sln_tree* tree, struct sln_node* leaf) 41 { 42 size_t vertices_range[2]; 43 res_T res = RES_OK; 44 45 /* Currently assume that we have only one line per leaf */ 46 ASSERT(leaf->range[0] == leaf->range[1]); 47 48 /* Line meshing */ 49 res = line_mesh(tree, leaf->range[0], tree->args.nvertices_hint, vertices_range); 50 if(res != RES_OK) goto error; 51 52 /* Decimate the line mesh */ 53 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 54 vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type); 55 if(res != RES_OK) goto error; 56 57 /* Shrink the size of the vertices */ 58 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 59 60 /* Setup the leaf polyline */ 61 leaf->ivertex = vertices_range[0]; 62 leaf->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 63 64 exit: 65 return res; 66 error: 67 goto exit; 68 } 69 70 static INLINE double 71 eval_ka 72 (const struct sln_tree* tree, 73 const struct sln_node* node, 74 const double wavenumber) 75 { 76 struct sln_mesh mesh = SLN_MESH_NULL; 77 double ka = 0; 78 ASSERT(tree && node); 79 80 /* Whether the mesh to be constructed corresponds to the spectrum or its upper 81 * limit, use the node mesh to calculate the value of ka at a given wave 82 * number. Calculating the value from the node lines would take far too long*/ 83 SLN(node_get_mesh(tree, node, &mesh)); 84 ka = sln_mesh_eval(&mesh, wavenumber); 85 86 return ka; 87 } 88 89 static res_T 90 merge_children_polylines 91 (struct sln_tree* tree, 92 const size_t inode) 93 { 94 /* Helper constant */ 95 static const size_t NO_MORE_VERTEX = SIZE_MAX; 96 97 /* Polyline vertices */ 98 size_t children_ivtx[SLN_TREE_ARITY_MAX] = {0}; 99 size_t vertices_range[2] = {0,0}; 100 struct sln_vertex* vertices = NULL; 101 size_t ivtx = 0; 102 size_t nvertices = 0; 103 104 /* Miscellaneous */ 105 size_t nchildren = 0; 106 unsigned i = 0; 107 res_T res = RES_OK; 108 109 /* Pre-conditions */ 110 ASSERT(tree && inode < darray_node_size_get(&tree->nodes)); 111 ASSERT(tree->args.arity >= 2); 112 ASSERT(tree->args.arity <= SLN_TREE_ARITY_MAX); 113 114 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 115 116 nchildren = node_child_count(tree, inode, NULL); 117 118 /* Compute the number of vertices to be merged, 119 * i.e., the sum of vertices of the children */ 120 nvertices = 0; 121 FOR_EACH(i, 0, nchildren) { 122 const size_t ichild = inode + NODE(inode)->offset + i; 123 nvertices += NODE(ichild)->nvertices; 124 } 125 126 /* Define the vertices range of the merged polyline */ 127 vertices_range[0] = darray_vertex_size_get(&tree->vertices); 128 vertices_range[1] = vertices_range[0] + nvertices - 1/*inclusive bound*/; 129 130 /* Allocate the memory space to store the new polyline */ 131 res = darray_vertex_resize(&tree->vertices, vertices_range[1]+1); 132 if(res != RES_OK) { 133 ERROR(tree->sln, "Error in merging polylines -- %s.\n", res_to_cstr(res)); 134 goto error; 135 } 136 vertices = darray_vertex_data_get(&tree->vertices); 137 138 /* Initialize the vertex index list. For each child, the initial value 139 * corresponds to the index of its first vertex. This index will be 140 * incremented as vertices are merged into the parent polyline. */ 141 FOR_EACH(i, 0, nchildren) { 142 const size_t ichild = inode + NODE(inode)->offset + i; 143 children_ivtx[i] = NODE(ichild)->ivertex; 144 } 145 146 FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1/*inclusive bound*/) { 147 double ka = 0; 148 double nu = INF; 149 150 /* The number of vertices corresponding to the current wave number for which 151 * the parent ka is calculated. It is at least equal to one, since this nu 152 * is defined by the child vertices, but may be greater if multiple children 153 * share the same vertex, i.e., a ka value calculated for the same nu */ 154 unsigned nvertices_merged = 0; 155 156 /* Find the minimum wave number among the vertices of the child vertices 157 * that are candidates for merging */ 158 FOR_EACH(i, 0, nchildren) { 159 const size_t child_ivtx = children_ivtx[i]; 160 if(child_ivtx != NO_MORE_VERTEX) { 161 nu = MMIN(nu, vertices[child_ivtx].wavenumber); 162 } 163 } 164 ASSERT(nu != INF); /* At least one vertex must have been found */ 165 166 /* Compute the value of ka at the wave number determined above */ 167 FOR_EACH(i, 0, nchildren) { 168 const size_t child_ivtx = children_ivtx[i]; 169 const struct sln_node* child = NODE(inode + NODE(inode)->offset + i); 170 171 if(child_ivtx == NO_MORE_VERTEX 172 || nu != vertices[child_ivtx].wavenumber) { 173 /* The wave number does not correspond to a vertex in the current 174 * child's mesh. Therefore, its contribution to the parent node's ka is 175 * computed */ 176 ka += eval_ka(tree, child, nu); 177 178 } else { 179 /* The wave number is the one for which the child node stores a ka 180 * value. Add it to the parent node's ka value and designate the child's 181 * next vertex as a candidate for merging into the parent. The exception 182 * is when all vertices of the child have already been merged. In this 183 * case, report that the child no longer has any candidate vertices */ 184 ka += vertices[child_ivtx].ka; 185 ++children_ivtx[i]; 186 if(children_ivtx[i] >= child->ivertex + child->nvertices) { 187 children_ivtx[i] = NO_MORE_VERTEX; 188 } 189 ++nvertices_merged; /* Record that a vertex has been merged */ 190 } 191 } 192 193 /* Setup the parent vertex */ 194 vertices[ivtx].wavenumber = (float)nu; 195 vertices[ivtx].ka = (float)ka; 196 197 /* If multiple child vertices have been merged, then a single wave number 198 * corresponds to a vertex with multiple children. The number of parent 199 * vertices is therefore no longer the sum of the number of its children's 200 * vertices, since some vertices are duplicated. Hence the following 201 * adjustment, which removes the duplicate vertices. */ 202 vertices_range[1] -= (nvertices_merged-1); 203 } 204 205 /* Decimate the resulting polyline */ 206 res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices), 207 vertices_range, (float)tree->args.mesh_decimation_err, tree->args.mesh_type); 208 if(res != RES_OK) goto error; 209 210 /* Shrink the size of the vertices */ 211 darray_vertex_resize(&tree->vertices, vertices_range[1] + 1); 212 213 /* Setup the node polyline */ 214 NODE(inode)->ivertex = vertices_range[0]; 215 NODE(inode)->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1); 216 217 /* It is necessary to ensure that the recorded vertices define a polyline 218 * along which any value (calculated by linear interpolation) is well above 219 * the sum of the corresponding values of the polylines to be merged. However, 220 * although this is guaranteed by definition for the vertices of the polyline, 221 * numerical uncertainty may nevertheless introduce errors that violate this 222 * criterion. Hence the following adjustment, which slightly increases the ka 223 * of the mesh so as to guarantee this constraint between the mesh of a node 224 * and that of its children */ 225 if(tree->args.mesh_type == SLN_MESH_UPPER) { 226 FOR_EACH(ivtx, vertices_range[0], vertices_range[1]+1/*inclusive bound*/) { 227 const double ka = vertices[ivtx].ka; 228 vertices[ivtx].ka = (float)(ka + MMAX(ka*1e-2, 1e-6)); 229 } 230 } 231 232 #undef NODE 233 234 exit: 235 return res; 236 error: 237 goto exit; 238 } 239 240 static res_T 241 build_polylines(struct sln_tree* tree) 242 { 243 /* Stack */ 244 #define STACK_SIZE (SLN_TREE_DEPTH_MAX*SLN_TREE_ARITY_MAX) 245 size_t stack[STACK_SIZE]; 246 size_t istack = 0; 247 248 /* Progress */ 249 size_t nnodes_total = 0; 250 size_t nnodes_processed = 0; 251 int progress = 0; 252 253 /* Miscellaneous */ 254 size_t inode = 0; 255 res_T res = RES_OK; 256 257 ASSERT(tree); 258 259 #define LOG_MSG "Meshing: %3d%%\n" 260 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 261 #define IS_LEAF(Id) (NODE(Id)->offset == 0) 262 263 nnodes_total = darray_node_size_get(&tree->nodes); 264 265 /* Print that nothing has been done yet */ 266 INFO(tree->sln, LOG_MSG, progress); 267 268 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 269 stack[istack++] = SIZE_MAX; 270 271 inode = 0; /* Root node */ 272 while(inode != SIZE_MAX) { 273 const size_t istack_saved = istack; 274 275 if(IS_LEAF(inode)) { 276 res = build_leaf_polyline(tree, NODE(inode)); 277 if(res != RES_OK) goto error; 278 279 inode = stack[--istack]; /* Pop the next node */ 280 281 } else { 282 const size_t ichild0 = inode + NODE(inode)->offset + 0; 283 const size_t nchildren = node_child_count(tree, inode, NULL); 284 size_t i = 0; 285 286 /* Child nodes have their polyline created */ 287 if(NODE(ichild0)->nvertices) { 288 #ifndef NDEBUG 289 /* Check that all children have their polylines created */ 290 FOR_EACH(i, 1, nchildren) { 291 const size_t ichild = inode + NODE(inode)->offset + i; 292 ASSERT(NODE(ichild)->nvertices != 0); 293 } 294 #endif 295 res = merge_children_polylines(tree, inode); 296 if(res != RES_OK) goto error; 297 inode = stack[--istack]; /* Pop the next node */ 298 299 /* Child nodes have NOT their polyline created */ 300 } else { 301 ASSERT(istack < STACK_SIZE - (nchildren - 1/*ichild0*/ + 1/*inode*/)); 302 stack[istack++] = inode; /* Push the current node */ 303 304 /* Push the children except the 1st one */ 305 FOR_EACH(i, 1, nchildren) { 306 const size_t ichild = inode + NODE(inode)->offset + i; 307 ASSERT(NODE(ichild)->nvertices == 0); 308 stack[istack++] = ichild; 309 } 310 311 /* Recursively build the polyline of the 1st child */ 312 inode = ichild0; 313 } 314 } 315 316 /* Handle progression bar */ 317 if(istack < istack_saved) { 318 int pcent = 0; 319 320 nnodes_processed += istack_saved - istack; 321 ASSERT(nnodes_processed <= nnodes_total); 322 323 /* Print progress update */ 324 pcent = (int)((double)nnodes_processed*100.0/(double)nnodes_total+0.5); 325 if(pcent/10 > progress/10) { 326 progress = pcent; 327 INFO(tree->sln, LOG_MSG, progress); 328 } 329 } 330 } 331 ASSERT(nnodes_processed == nnodes_total); 332 333 #undef NODE 334 #undef IS_LEAF 335 #undef LOG_MSG 336 #undef STACK_SIZE 337 338 exit: 339 return res; 340 error: 341 goto exit; 342 } 343 344 static res_T 345 partition_lines(struct sln_tree* tree) 346 { 347 /* Stack */ 348 #define STACK_SIZE (SLN_TREE_DEPTH_MAX*(SLN_TREE_ARITY_MAX-1)) 349 size_t stack[STACK_SIZE]; 350 size_t istack = 0; 351 352 /* Progress */ 353 size_t nlines = 0; 354 size_t nlines_total = 0; 355 size_t nlines_processed = 0; 356 int progress = 0; 357 358 /* Miscellaneous */ 359 size_t inode = 0; 360 res_T res = RES_OK; 361 362 ASSERT(tree); /* Pre-condition */ 363 364 SHTR(line_list_get_size(tree->args.lines, &nlines)); 365 366 nlines_total = nlines; 367 nlines_processed = 0; 368 369 #define LOG_MSG "Partitioning: %3d%%\n" 370 #define NODE(Id) (darray_node_data_get(&tree->nodes) + (Id)) 371 #define CREATE_NODE { \ 372 res = darray_node_push_back(&tree->nodes, &SLN_NODE_NULL); \ 373 if(res != RES_OK) goto error; \ 374 } (void)0 375 376 CREATE_NODE; /* Alloc the root node */ 377 378 /* Setup the root node */ 379 NODE(0)->range[0] = 0; 380 NODE(0)->range[1] = nlines - 1; 381 382 /* Push back SIZE_MAX which, once pop up, will mark the end of recursion */ 383 stack[istack++] = SIZE_MAX; 384 385 /* Print that although nothing has been done yet, 386 * the calculation is nevertheless in progress */ 387 INFO(tree->sln, LOG_MSG, progress); 388 389 inode = 0; /* Root node */ 390 while(inode != SIZE_MAX) { 391 /* #lines into the node */ 392 nlines = NODE(inode)->range[1] - NODE(inode)->range[0] + 1; 393 394 /* Make a leaf */ 395 if(nlines <= MAX_NLINES_PER_LEAF) { 396 int pcent = 0; 397 398 NODE(inode)->offset = 0; 399 inode = stack[--istack]; /* Pop the next node */ 400 401 ASSERT(nlines_processed + nlines <= nlines_total); 402 nlines_processed += nlines; 403 404 /* Print progress update */ 405 pcent = (int)((double)nlines_processed*100.0/(double)nlines_total+0.5); 406 if(pcent/10 > progress/10) { 407 progress = pcent; 408 INFO(tree->sln, LOG_MSG, progress); 409 } 410 411 /* Split the node */ 412 } else { 413 size_t node_range[2] = {0,0}; 414 size_t nlines_per_child = 0; 415 size_t nchildren = 0; 416 size_t i = 0; 417 418 /* Compute how the lines of a node are distributed among its children, 419 * as well as the number of children that node has */ 420 nchildren = node_child_count(tree, inode, &nlines_per_child); 421 422 /* Calculate the index of the first child */ 423 size_t ichildren = darray_node_size_get(&tree->nodes); 424 425 ASSERT(nchildren <= tree->args.arity); 426 ASSERT(ichildren > inode); 427 428 node_range[0] = NODE(inode)->range[0]; 429 node_range[1] = NODE(inode)->range[1]; 430 431 /* Define the offset from the current node to its children */ 432 NODE(inode)->offset = ui64_to_ui32((uint64_t)(ichildren - inode)); 433 434 FOR_EACH(i, 0, nchildren) { 435 CREATE_NODE; 436 437 /* Set the range of lines line for the newly created child. Note that 438 * the boundaries of the range are inclusive, which is why 1 is 439 * subtracted to the upper bound */ 440 NODE(ichildren+i)->range[0] = node_range[0] + nlines_per_child*(i+0); 441 NODE(ichildren+i)->range[1] = node_range[0] + nlines_per_child*(i+1)-1; 442 443 /* Check that the child's lines are a subset of the parent's lines. Note 444 * that the upper bound of the child's line interval is not checked, as 445 * it may temporarily exceed the parent's interval since, up to this 446 * point, it is assumed that all children contain the same number of 447 * lines. This is no longer true if the number of lines in the parent is 448 * not a multiple of the computed number of lines per child. The upper 449 * bound of the last child's range is therefore readjusted at the end of 450 * the loop (see below) */ 451 ASSERT(NODE(ichildren+i)->range[0] <= node_range[1]); 452 } 453 /* Readjust the upper bound of the last child to ensure that its line 454 * interval is fully contained within its parent node. See above for an 455 * explanation of why */ 456 NODE(ichildren + nchildren - 1)->range[1] = node_range[1]; 457 458 inode = ichildren; /* Make the first child the current node */ 459 460 ASSERT(istack < STACK_SIZE - 1/*1st child*/); 461 FOR_EACH(i, 1, nchildren) { 462 stack[istack++] = ichildren + i; /* Push the other children */ 463 } 464 } 465 } 466 ASSERT(nlines_processed == nlines_total); 467 468 #undef NODE 469 #undef CREATE_NODE 470 #undef STACK_SIZE 471 exit: 472 return res; 473 error: 474 darray_node_clear(&tree->nodes); 475 goto exit; 476 } 477 478 /******************************************************************************* 479 * Local functions 480 ******************************************************************************/ 481 res_T 482 tree_build(struct sln_tree* tree) 483 { 484 res_T res = RES_OK; 485 486 if((res = partition_lines(tree)) != RES_OK) goto error; 487 if((res = build_polylines(tree)) != RES_OK) goto error; 488 489 exit: 490 return res; 491 error: 492 goto exit; 493 }