sck.c (30933B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 CNRS 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200809L /* mmap support */ 18 #define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */ 19 #define _BSD_SOURCE 1 /* MAP_POPULATE for glibc < 2.19 */ 20 21 #include "sck.h" 22 #include "sck_c.h" 23 #include "sck_log.h" 24 25 #include <rsys/algorithm.h> 26 #include <rsys/cstr.h> 27 28 #include <unistd.h> /* sysconf support */ 29 30 #include <errno.h> 31 #include <string.h> 32 #include <sys/mman.h> /* mmap */ 33 #include <sys/stat.h> /* fstat */ 34 35 /******************************************************************************* 36 * Helper functions 37 ******************************************************************************/ 38 static INLINE int 39 is_stdin(FILE* stream) 40 { 41 struct stat stream_buf; 42 struct stat stdin_buf; 43 ASSERT(stream); 44 CHK(fstat(fileno(stream), &stream_buf) == 0); 45 CHK(fstat(STDIN_FILENO, &stdin_buf) == 0); 46 return stream_buf.st_dev == stdin_buf.st_dev; 47 } 48 49 static INLINE res_T 50 check_sck_create_args(const struct sck_create_args* args) 51 { 52 /* Nothing to check. Only return RES_BAD_ARG if args is NULL */ 53 return args ? RES_OK : RES_BAD_ARG; 54 } 55 56 static INLINE res_T 57 check_sck_load_args(const struct sck_load_args* args) 58 { 59 if(!args || !args->path) return RES_BAD_ARG; 60 return RES_OK; 61 } 62 63 static INLINE res_T 64 check_sck_load_stream_args 65 (struct sck* sck, 66 const struct sck_load_stream_args* args) 67 { 68 if(!args || !args->stream || !args->name) return RES_BAD_ARG; 69 if(args->memory_mapping && is_stdin(args->stream)) { 70 log_err(sck, "%s: unable to use memory mapping on data loaded from stdin\n", 71 args->name); 72 return RES_BAD_ARG; 73 } 74 return RES_OK; 75 } 76 77 static void 78 reset_sck(struct sck* sck) 79 { 80 ASSERT(sck); 81 sck->pagesize = 0; 82 sck->nnodes = 0; 83 darray_band_purge(&sck->bands); 84 str_clear(&sck->name); 85 } 86 87 static INLINE int 88 cmp_dbl(const void* a, const void* b) 89 { 90 const double d0 = *((const double*)a); 91 const double d1 = *((const double*)b); 92 return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0); 93 } 94 95 static res_T 96 read_quad_pt 97 (struct sck* sck, 98 struct quad_pt* quad_pt, 99 FILE* stream, 100 size_t iband, 101 size_t iquad_pt) 102 { 103 res_T res = RES_OK; 104 ASSERT(sck && quad_pt); 105 106 quad_pt->band = darray_band_data_get(&sck->bands) + iband; 107 108 if(fread(&quad_pt->weight, sizeof(quad_pt->weight), 1, stream) != 1) { 109 log_err(sck, 110 "%s: band %lu: quadrature point %lu: could not read the weight.\n", 111 sck_get_name(sck), iband, iquad_pt); 112 res = RES_IO_ERR; 113 goto error; 114 } 115 116 if(quad_pt->weight < 0) { 117 log_err(sck, 118 "%s: band %lu: quadrature point %lu: invalid weight %g.\n", 119 sck_get_name(sck), iband, iquad_pt, quad_pt->weight); 120 res = RES_BAD_ARG; 121 goto error; 122 123 } 124 125 quad_pt->ka_list = NULL; 126 127 exit: 128 return res; 129 error: 130 goto exit; 131 } 132 133 static res_T 134 read_band 135 (struct sck* sck, 136 struct band* band, 137 FILE* stream) 138 { 139 double cumul; 140 size_t iquad_pt; 141 size_t iband; 142 uint64_t nquad_pts; 143 res_T res = RES_OK; 144 ASSERT(sck && band); 145 146 band->sck = sck; 147 iband = (size_t)(band - darray_band_cdata_get(&sck->bands)); 148 149 /* Read band definition */ 150 #define READ(Var, Name) { \ 151 if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ 152 log_err(sck, "%s: band %lu: could not read the %s.\n", \ 153 sck_get_name(sck), (unsigned long)iband, (Name)); \ 154 res = RES_IO_ERR; \ 155 goto error; \ 156 } \ 157 } (void)0 158 READ(&band->low, "band lower bound"); 159 READ(&band->upp, "band upper bound"); 160 READ(&nquad_pts, "#quadrature points"); 161 #undef READ 162 163 /* Check band description */ 164 if(band->low < 0 || band->low > band->upp) { 165 log_err(sck, 166 "%s: band %lu: invalid band range [%g, %g].\n", 167 sck_get_name(sck), (unsigned long)iband, band->low, band->upp); 168 res = RES_BAD_ARG; 169 goto error; 170 } 171 if(nquad_pts == 0) { 172 log_err(sck, 173 "%s: band %lu: invalid number fo quadrature points (#points=%lu).\n", 174 sck_get_name(sck), (unsigned long)iband, (unsigned long)nquad_pts); 175 res = RES_BAD_ARG; 176 goto error; 177 } 178 179 /* Allocate the list of quadrature points */ 180 res = darray_quad_pt_resize(&band->quad_pts, nquad_pts); 181 if(res != RES_OK) { 182 log_err(sck, 183 "%s: band %lu: could not allocate the list of quadrature points " 184 "(#points=%lu).\n", 185 sck_get_name(sck), (unsigned long)iband, (unsigned long)nquad_pts); 186 goto error; 187 } 188 189 /* Allocate the cumulative used to sample the quadrature points */ 190 res = darray_double_resize(&band->quad_pts_cumul, nquad_pts); 191 if(res != RES_OK) { 192 log_err(sck, 193 "%s: band %lu: could not allocate the cumulative of quadrature points " 194 "(#points=%lu).\n", 195 sck_get_name(sck), (unsigned long)iband, (unsigned long)nquad_pts); 196 goto error; 197 } 198 199 /* Read the data of the quadrature points of the band */ 200 cumul = 0; 201 FOR_EACH(iquad_pt, 0, nquad_pts) { 202 struct quad_pt* quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; 203 204 /* Read current quadrature point */ 205 res = read_quad_pt(sck, quad_pt, stream, iband, iquad_pt); 206 if(res != RES_OK) goto error; 207 208 /* Compute the quadrature point cumulative */ 209 cumul += quad_pt->weight; 210 darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] = cumul; 211 } 212 213 /* The weights are not normalized */ 214 if(!eq_eps(cumul, 1.0, 1.e-6)) { 215 log_warn(sck, 216 "%s: band %lu: the weights of the quadrature points are not normalised.\n", 217 sck_get_name(sck), (unsigned long)iband); 218 219 /* Renormalize the cumulative */ 220 FOR_EACH(iquad_pt, 0, nquad_pts) { 221 darray_double_data_get(&band->quad_pts_cumul)[iquad_pt] /= cumul; 222 } 223 } 224 225 /* Handle imprecision issue */ 226 darray_double_data_get(&band->quad_pts_cumul)[nquad_pts-1] = 1.0; 227 228 exit: 229 return res; 230 error: 231 goto exit; 232 } 233 234 static res_T 235 map_data 236 (struct sck* sck, 237 const int fd, /* File descriptor */ 238 const size_t filesz, /* Overall filesize */ 239 const char* data_name, 240 const off_t offset, /* Offset of the data into file */ 241 const size_t map_len, 242 void** out_map) /* Lenght of the data to map */ 243 { 244 void* map = NULL; 245 res_T res = RES_OK; 246 ASSERT(sck && filesz && data_name && map_len && out_map); 247 ASSERT(IS_ALIGNED((size_t)offset, (size_t)sck->pagesize)); 248 249 if((size_t)offset + map_len > filesz) { 250 log_err(sck, "%s: the %s to map exceed the file size\n", 251 sck_get_name(sck), data_name); 252 res = RES_IO_ERR; 253 goto error; 254 } 255 256 map = mmap(NULL, map_len, PROT_READ, MAP_PRIVATE|MAP_POPULATE, fd, offset); 257 if(map == MAP_FAILED) { 258 log_err(sck, "%s: could not map the %s -- %s\n", 259 sck_get_name(sck), data_name, strerror(errno)); 260 res = RES_IO_ERR; 261 goto error; 262 } 263 264 exit: 265 *out_map = map; 266 return res; 267 error: 268 if(map == MAP_FAILED) map = NULL; 269 goto exit; 270 } 271 272 static res_T 273 map_file(struct sck* sck, FILE* stream) 274 { 275 size_t filesz; 276 size_t map_len; 277 size_t iband; 278 off_t offset; 279 res_T res = RES_OK; 280 ASSERT(sck && stream); 281 282 /* Compute the length in bytes of the k to map for each band/quadrature point */ 283 map_len = ALIGN_SIZE(sck->nnodes * sizeof(float), sck->pagesize); 284 285 /* Compute the offset toward the 1st list of radiative coefficients */ 286 offset = ftell(stream); 287 offset = (off_t)ALIGN_SIZE((uint64_t)offset, sck->pagesize); 288 289 /* Retrieve the overall filesize */ 290 fseek(stream, 0, SEEK_END); 291 filesz = (size_t)ftell(stream); 292 293 FOR_EACH(iband, 0, darray_band_size_get(&sck->bands)) { 294 struct band* band = NULL; 295 size_t iquad_pt; 296 297 band = darray_band_data_get(&sck->bands) + iband; 298 band->map_len = map_len; 299 300 /* Mapping per band scattering coefficients */ 301 res = map_data(sck, fileno(stream), filesz, "scattering coefficients", 302 offset, band->map_len, (void**)&band->ks_list); 303 if(res != RES_OK) { 304 log_err(sck, 305 "%s: data mapping error for band %lu\n", 306 sck_get_name(sck), (unsigned long)iband); 307 goto error; 308 } 309 310 offset = (off_t)((size_t)offset + map_len); 311 312 FOR_EACH(iquad_pt, 0, darray_quad_pt_size_get(&band->quad_pts)) { 313 struct quad_pt* quad_pt = NULL; 314 315 /* Mapping absorption coefficients per band and quadrature point */ 316 quad_pt = darray_quad_pt_data_get(&band->quad_pts)+iquad_pt; 317 quad_pt->map_len = map_len; 318 319 res = map_data(sck, fileno(stream), filesz, "correlated Ka", 320 offset, quad_pt->map_len, (void**)&quad_pt->ka_list); 321 if(res != RES_OK) { 322 log_err(sck, 323 "%s: data mapping error for band %lu quadrature point %lu\n", 324 sck_get_name(sck), (unsigned long)iband, (unsigned long)iquad_pt); 325 goto error; 326 } 327 328 offset = (off_t)((size_t)offset + map_len); 329 } 330 } 331 332 exit: 333 return res; 334 error: 335 goto exit; 336 } 337 338 static res_T 339 read_padding(FILE* stream, const size_t padding) 340 { 341 char chunk[1024]; 342 size_t remaining_nbytes = padding; 343 344 while(remaining_nbytes) { 345 const size_t nbytes = MMIN(sizeof(chunk), remaining_nbytes); 346 if(fread(chunk, 1, nbytes, stream) != nbytes) return RES_IO_ERR; 347 remaining_nbytes -= nbytes; 348 } 349 return RES_OK; 350 } 351 352 /* Return the size in bytes of the data layout and band descriptors */ 353 static size_t 354 compute_sizeof_header(struct sck* sck) 355 { 356 size_t sizeof_header = 0; 357 size_t iband = 0; 358 size_t nbands = 0; 359 ASSERT(sck); 360 361 sizeof_header = 362 sizeof(uint64_t) /* pagesize */ 363 + sizeof(uint64_t) /* #bands */ 364 + sizeof(uint64_t);/* #nodes */ 365 366 nbands = darray_band_size_get(&sck->bands); 367 FOR_EACH(iband, 0, nbands) { 368 const struct band* band = darray_band_cdata_get(&sck->bands)+iband; 369 const size_t nquad_pts = darray_quad_pt_size_get(&band->quad_pts); 370 const size_t sizeof_band_desc = 371 sizeof(double) /* Lower bound */ 372 + sizeof(double) /* Upper bound */ 373 + sizeof(uint64_t) /* #quad_pts */ 374 + sizeof(double) * nquad_pts; /* Weights */ 375 376 sizeof_header += sizeof_band_desc; 377 } 378 return sizeof_header; 379 } 380 381 static res_T 382 load_data 383 (struct sck* sck, 384 FILE* stream, 385 const char* data_name, 386 float** out_data) 387 { 388 float* data = NULL; 389 res_T res = RES_OK; 390 ASSERT(sck && stream && data_name && out_data); 391 392 data = MEM_ALLOC(sck->allocator, sizeof(float)*sck->nnodes); 393 if(!data) { 394 res = RES_MEM_ERR; 395 log_err(sck, "%s: could not allocate the %s -- %s\n", 396 sck_get_name(sck), data_name, res_to_cstr(res)); 397 goto error; 398 } 399 400 if(fread(data, sizeof(float), sck->nnodes, stream) != sck->nnodes) { 401 res = RES_IO_ERR; 402 log_err(sck, "%s: could not read the %s -- %s\n", 403 sck_get_name(sck), data_name, res_to_cstr(res)); 404 goto error; 405 } 406 407 exit: 408 *out_data = data; 409 return res; 410 error: 411 if(data) { MEM_RM(sck->allocator, data); data = NULL; } 412 goto exit; 413 } 414 415 static res_T 416 load_file(struct sck* sck, FILE* stream) 417 { 418 size_t padding_bytes; 419 size_t sizeof_header; 420 size_t sizeof_k_list; 421 size_t iband; 422 size_t nbands; 423 res_T res = RES_OK; 424 425 sizeof_header = compute_sizeof_header(sck); 426 sizeof_k_list = sizeof(float)*sck->nnodes; 427 428 padding_bytes = ALIGN_SIZE(sizeof_header, sck->pagesize) - sizeof_header; 429 if((res = read_padding(stream, padding_bytes)) != RES_OK) goto error; 430 431 /* Calculate the padding between the lists of radiative coefficients. Note 432 * that this padding is the same between each list */ 433 padding_bytes = ALIGN_SIZE(sizeof_k_list, sck->pagesize) - sizeof_k_list; 434 435 nbands = darray_band_size_get(&sck->bands); 436 FOR_EACH(iband, 0, nbands) { 437 struct band* band = NULL; 438 size_t iquad_pt = 0; 439 440 band = darray_band_data_get(&sck->bands) + iband; 441 ASSERT(!band->ks_list && band->sck == sck); 442 443 /* Loading per band scattering coefficients */ 444 res = load_data(sck, stream, "scattering coefficients", &band->ks_list); 445 if(res != RES_OK) { 446 log_err(sck, 447 "%s: data loading error for band %lu\n", 448 sck_get_name(sck), (unsigned long)iband); 449 goto error; 450 } 451 452 if((res = read_padding(stream, padding_bytes)) != RES_OK) goto error; 453 454 FOR_EACH(iquad_pt, 0, darray_quad_pt_size_get(&band->quad_pts)) { 455 struct quad_pt* quad_pt = NULL; 456 457 quad_pt = darray_quad_pt_data_get(&band->quad_pts) + iquad_pt; 458 459 /* Loading absorption coefficients per band and quadrature point */ 460 res = load_data(sck, stream, "correlated Ka", &quad_pt->ka_list); 461 if(res != RES_OK) { 462 log_err(sck, 463 "%s: data loading error for band %lu quadrature point %lu\n", 464 sck_get_name(sck), (unsigned long)iband, (unsigned long)iquad_pt); 465 goto error; 466 } 467 468 if((res = read_padding(stream, padding_bytes)) != RES_OK) goto error; 469 } 470 } 471 472 exit: 473 return res; 474 error: 475 goto exit; 476 } 477 478 static res_T 479 load_stream(struct sck* sck, const struct sck_load_stream_args* args) 480 { 481 size_t iband; 482 uint64_t nbands; 483 res_T res = RES_OK; 484 ASSERT(sck && check_sck_load_stream_args(sck, args) == RES_OK); 485 486 reset_sck(sck); 487 488 res = str_set(&sck->name, args->name); 489 if(res != RES_OK) { 490 log_err(sck, "%s: unable to duplicate path to loaded data or stream name\n", 491 args->name); 492 goto error; 493 } 494 495 /* Read file header */ 496 if(fread(&sck->pagesize, sizeof(sck->pagesize), 1, args->stream) != 1) { 497 if(ferror(args->stream)) { 498 log_err(sck, "%s: could not read the pagesize.\n", sck_get_name(sck)); 499 } 500 res = RES_IO_ERR; 501 goto error; 502 } 503 #define READ(Var, Name) { \ 504 if(fread((Var), sizeof(*(Var)), 1, args->stream) != 1) { \ 505 log_err(sck, "%s: could not read the %s.\n", sck_get_name(sck), (Name)); \ 506 res = RES_IO_ERR; \ 507 goto error; \ 508 } \ 509 } (void)0 510 READ(&nbands, "number of bands"); 511 READ(&sck->nnodes, "number of nodes"); 512 #undef READ 513 514 /* Check band description */ 515 if(!IS_ALIGNED(sck->pagesize, sck->pagesize_os)) { 516 log_err(sck, 517 "%s: invalid page size %lu. The page size attribute must be aligned on " 518 "the page size of the operating system (%lu).\n", 519 sck_get_name(sck), sck->pagesize, (unsigned long)sck->pagesize_os); 520 res = RES_BAD_ARG; 521 goto error; 522 } 523 524 if(!nbands) { 525 log_err(sck, "%s: invalid number of bands %lu.\n", 526 sck_get_name(sck), (unsigned long)nbands); 527 res = RES_BAD_ARG; 528 goto error; 529 } 530 if(!sck->nnodes) { 531 log_err(sck, "%s: invalid number of nodes %lu.\n", 532 sck_get_name(sck), (unsigned long)sck->nnodes); 533 res = RES_BAD_ARG; 534 goto error; 535 } 536 537 /* Allocate the bands */ 538 res = darray_band_resize(&sck->bands, nbands); 539 if(res != RES_OK) { 540 log_err(sck, "%s: could not allocate the list of bands (#bands=%lu).\n", 541 sck_get_name(sck), (unsigned long)nbands); 542 goto error; 543 } 544 545 /* Read the band description */ 546 FOR_EACH(iband, 0, nbands) { 547 struct band* band = darray_band_data_get(&sck->bands) + iband; 548 res = read_band(sck, band, args->stream); 549 if(res != RES_OK) goto error; 550 if(iband > 0 && band[0].low < band[-1].upp) { 551 log_err(sck, 552 "%s: bands must be sorted in ascending order and must not " 553 "overlap (band %lu in [%g, %g[ nm; band %lu in [%g, %g[ nm).\n", 554 sck_get_name(sck), 555 (unsigned long)(iband-1), band[-1].low, band[-1].upp, 556 (unsigned long)(iband), band[ 0].low, band[ 0].upp); 557 res = RES_BAD_ARG; 558 goto error; 559 } 560 } 561 562 if(args->memory_mapping) { 563 res = map_file(sck, args->stream); 564 if(res != RES_OK) goto error; 565 } else { 566 res = load_file(sck, args->stream); 567 if(res != RES_OK) goto error; 568 } 569 570 exit: 571 return res; 572 error: 573 reset_sck(sck); 574 goto exit; 575 } 576 577 static INLINE int 578 cmp_band(const void* key, const void* item) 579 { 580 const struct band* band = item; 581 double wnum; 582 ASSERT(key && item); 583 wnum = *(double*)key; 584 585 if(wnum < band->low) { 586 return -1; 587 } else if(wnum >= band->upp) { 588 return +1; 589 } else { 590 return 0; 591 } 592 } 593 594 static INLINE void 595 hash_quad_pt 596 (struct sha256_ctx* ctx, 597 const struct sck_quad_pt* pt, 598 const size_t nnodes) 599 { 600 ASSERT(ctx && pt); 601 #define HASH(Var, Nb) sha256_ctx_update(ctx, (const char*)(Var), sizeof(*Var)*(Nb)) 602 HASH(&pt->weight, 1); 603 HASH(&pt->id, 1); 604 HASH(pt->ka_list, nnodes); 605 #undef HASH 606 } 607 608 static INLINE void 609 hash_band(struct sha256_ctx* ctx, const struct sck_band* band) 610 { 611 size_t iquad_pt; 612 ASSERT(ctx && band); 613 614 #define HASH(Var, Nb) sha256_ctx_update(ctx, (const char*)(Var), sizeof(*Var)*(Nb)) 615 HASH(&band->lower, 1); 616 HASH(&band->upper, 1); 617 HASH(&band->quad_pts_count, 1); 618 HASH(band->ks_list, band->sck__->nnodes); 619 #undef HASH 620 621 FOR_EACH(iquad_pt, 0, band->quad_pts_count) { 622 struct sck_quad_pt quad_pt; 623 SCK(band_get_quad_pt(band, iquad_pt, &quad_pt)); 624 hash_quad_pt(ctx, &quad_pt, band->sck__->nnodes); 625 } 626 } 627 628 static void 629 release_sck(ref_T* ref) 630 { 631 struct sck* sck; 632 ASSERT(ref); 633 sck = CONTAINER_OF(ref, struct sck, ref); 634 if(sck->logger == &sck->logger__) logger_release(&sck->logger__); 635 darray_band_release(&sck->bands); 636 str_release(&sck->name); 637 MEM_RM(sck->allocator, sck); 638 } 639 640 /******************************************************************************* 641 * Exported functions 642 ******************************************************************************/ 643 res_T 644 sck_create 645 (const struct sck_create_args* args, 646 struct sck** out_sck) 647 { 648 struct sck* sck = NULL; 649 struct mem_allocator* allocator = NULL; 650 res_T res = RES_OK; 651 652 if(!out_sck) { res = RES_BAD_ARG; goto error; } 653 res = check_sck_create_args(args); 654 if(res != RES_OK) goto error; 655 656 allocator = args->allocator ? args->allocator : &mem_default_allocator; 657 sck = MEM_CALLOC(allocator, 1, sizeof(*sck)); 658 if(!sck) { 659 if(args->verbose) { 660 #define ERR_STR "Could not allocate the Star-CK device.\n" 661 if(args->logger) { 662 logger_print(args->logger, LOG_ERROR, ERR_STR); 663 } else { 664 fprintf(stderr, MSG_ERROR_PREFIX ERR_STR); 665 } 666 #undef ERR_STR 667 } 668 res = RES_MEM_ERR; 669 goto error; 670 } 671 ref_init(&sck->ref); 672 sck->allocator = allocator; 673 sck->verbose = args->verbose; 674 sck->pagesize_os = (size_t)sysconf(_SC_PAGESIZE); 675 str_init(allocator, &sck->name); 676 darray_band_init(allocator, &sck->bands); 677 if(args->logger) { 678 sck->logger = args->logger; 679 } else { 680 setup_log_default(sck); 681 } 682 683 exit: 684 if(out_sck) *out_sck = sck; 685 return res; 686 error: 687 if(sck) { 688 SCK(ref_put(sck)); 689 sck = NULL; 690 } 691 goto exit; 692 } 693 694 res_T 695 sck_ref_get(struct sck* sck) 696 { 697 if(!sck) return RES_BAD_ARG; 698 ref_get(&sck->ref); 699 return RES_OK; 700 } 701 702 res_T 703 sck_ref_put(struct sck* sck) 704 { 705 if(!sck) return RES_BAD_ARG; 706 ref_put(&sck->ref, release_sck); 707 return RES_OK; 708 } 709 710 res_T 711 sck_load(struct sck* sck, const struct sck_load_args* args) 712 { 713 struct sck_load_stream_args stream_args = SCK_LOAD_STREAM_ARGS_NULL; 714 FILE* file = NULL; 715 res_T res = RES_OK; 716 717 if(!sck) { res = RES_BAD_ARG; goto error; } 718 res = check_sck_load_args(args); 719 if(res != RES_OK) goto error; 720 721 file = fopen(args->path, "r"); 722 if(!file) { 723 log_err(sck, "%s: error opening file `%s'.\n", FUNC_NAME, args->path); 724 res = RES_IO_ERR; 725 goto error; 726 } 727 728 stream_args.stream = file; 729 stream_args.name = args->path; 730 stream_args.memory_mapping = args->memory_mapping; 731 res = load_stream(sck, &stream_args); 732 if(res != RES_OK) goto error; 733 734 exit: 735 if(file) fclose(file); 736 return res; 737 error: 738 goto exit; 739 } 740 741 res_T 742 sck_load_stream(struct sck* sck, const struct sck_load_stream_args* args) 743 { 744 res_T res = RES_OK; 745 if(!sck) return RES_BAD_ARG; 746 res = check_sck_load_stream_args(sck, args); 747 if(res != RES_OK) return res; 748 return load_stream(sck, args); 749 } 750 751 res_T 752 sck_validate(const struct sck* sck) 753 { 754 size_t iband; 755 if(!sck) return RES_BAD_ARG; 756 757 FOR_EACH(iband, 0, sck_get_bands_count(sck)) { 758 struct sck_band band = SCK_BAND_NULL; 759 size_t inode; 760 size_t iquad_pt; 761 762 SCK(get_band(sck, iband, &band)); 763 764 /* Check band limits */ 765 if(band.lower != band.lower /* NaN? */ 766 || band.upper != band.upper) { /* NaN? */ 767 log_err(sck, "%s: invalid limits for band %lu: [%g, %g[\n", 768 sck_get_name(sck), (unsigned long)iband, band.lower, band.upper); 769 return RES_BAD_ARG; 770 } 771 772 /* Check scattering coefficients */ 773 FOR_EACH(inode, 0, sck_get_nodes_count(sck)) { 774 if(band.ks_list[inode] != band.ks_list[inode] /* NaN? */ 775 || band.ks_list[inode] < 0) { 776 log_err(sck, 777 "%s: invalid scattering coefficient for band %lu at node %lu: %g\n", 778 sck_get_name(sck), (unsigned long)iband, (unsigned long) inode, 779 band.ks_list[inode]); 780 return RES_BAD_ARG; 781 } 782 } 783 784 FOR_EACH(iquad_pt, 0, band.quad_pts_count) { 785 struct sck_quad_pt quad_pt = SCK_QUAD_PT_NULL; 786 787 SCK(band_get_quad_pt(&band, iquad_pt, &quad_pt)); 788 789 /* Check quadrature point weight */ 790 if(quad_pt.weight != quad_pt.weight /* NaN? */ 791 || quad_pt.weight <= 0) { 792 log_err(sck, 793 "%s: invalid weight for quadrature point %lu of band %lu: %g\n", 794 sck_get_name(sck), (unsigned long)iquad_pt, (unsigned long)iband, 795 quad_pt.weight); 796 return RES_BAD_ARG; 797 } 798 799 /* Check absorption coefficient */ 800 FOR_EACH(inode, 0, sck_get_nodes_count(sck)) { 801 if(quad_pt.ka_list[inode] != quad_pt.ka_list[inode] /* NaN? */ 802 || quad_pt.ka_list[inode] < 0) { 803 log_err(sck, 804 "%s: invalid absorption coefficient for quadrature point %lu " 805 "of band %lu at node %lu: %g\n", 806 sck_get_name(sck), (unsigned long)iquad_pt, (unsigned long)iband, 807 (unsigned long)inode, quad_pt.ka_list[inode]); 808 return RES_BAD_ARG; 809 } 810 } 811 } 812 } 813 return RES_OK; 814 } 815 816 size_t 817 sck_get_bands_count(const struct sck* sck) 818 { 819 ASSERT(sck); 820 return darray_band_size_get(&sck->bands); 821 } 822 823 size_t 824 sck_get_nodes_count(const struct sck* sck) 825 { 826 ASSERT(sck); 827 return sck->nnodes; 828 } 829 830 res_T 831 sck_get_band 832 (const struct sck* sck, 833 const size_t iband, 834 struct sck_band* sck_band) 835 { 836 const struct band* band = NULL; 837 res_T res = RES_OK; 838 839 if(!sck || !sck_band) { 840 res = RES_BAD_ARG; 841 goto error; 842 } 843 844 if(iband >= sck_get_bands_count(sck)) { 845 log_err(sck, "%s: invalid band index %lu.\n", 846 FUNC_NAME, (unsigned long)iband); 847 res = RES_BAD_ARG; 848 goto error; 849 } 850 851 band = darray_band_cdata_get(&sck->bands) + iband; 852 sck_band->lower = band->low; 853 sck_band->upper = band->upp; 854 sck_band->ks_list = band->ks_list; 855 sck_band->quad_pts_count = darray_quad_pt_size_get(&band->quad_pts); 856 sck_band->id = iband; 857 sck_band->sck__ = sck; 858 sck_band->band__ = band; 859 860 exit: 861 return res; 862 error: 863 goto exit; 864 } 865 866 res_T 867 sck_find_bands 868 (const struct sck* sck, 869 const double range[2], 870 size_t ibands[2]) 871 { 872 const struct band* bands = NULL; 873 const struct band* low = NULL; 874 const struct band* upp = NULL; 875 size_t nbands = 0; 876 res_T res = RES_OK; 877 878 if(!sck || !range || !ibands || range[0] > range[1]) { 879 res = RES_BAD_ARG; 880 goto error; 881 } 882 883 bands = darray_band_cdata_get(&sck->bands); 884 nbands = darray_band_size_get(&sck->bands); 885 886 low = search_lower_bound(range+0, bands, nbands, sizeof(*bands), cmp_band); 887 if(low) { 888 ibands[0] = (size_t)(low - bands); 889 } else { 890 /* The submitted range does not overlap any band */ 891 ibands[0] = SIZE_MAX; 892 ibands[1] = 0; 893 goto exit; 894 } 895 896 if(range[0] == range[1]) { /* No more to search */ 897 if(range[0] < low->low) { 898 /* The wavelength is not included in any band */ 899 ibands[0] = SIZE_MAX; 900 ibands[1] = 0; 901 } else { 902 ASSERT(range[0] < low->upp); 903 ibands[1] = ibands[0]; 904 } 905 goto exit; 906 } 907 908 upp = search_lower_bound(range+1, bands, nbands, sizeof(*bands), cmp_band); 909 910 /* The submitted range overlaps the remaining bands */ 911 if(!upp) { 912 ibands[1] = nbands - 1; 913 914 /* The upper band includes range[1] */ 915 } else if(upp->low <= range[1]) { 916 ibands[1] = (size_t)(upp - bands); 917 918 /* The upper band is greater than range[1] and therefre must be rejected */ 919 } else if(upp->low > range[1]) { 920 if(upp != bands) { 921 ibands[1] = (size_t)(upp - bands - 1); 922 } else { 923 ibands[0] = SIZE_MAX; 924 ibands[1] = 0; 925 } 926 } 927 928 exit: 929 return res; 930 error: 931 goto exit; 932 } 933 934 res_T 935 sck_band_get_quad_pt 936 (const struct sck_band* sck_band, 937 const size_t iquad_pt, 938 struct sck_quad_pt* sck_quad_pt) 939 { 940 const struct band* band = NULL; 941 const struct quad_pt* quad_pt = NULL; 942 res_T res = RES_OK; 943 944 if(!sck_band || !sck_quad_pt) { 945 res = RES_BAD_ARG; 946 goto error; 947 } 948 949 band = sck_band->band__; 950 if(iquad_pt >= sck_band->quad_pts_count) { 951 log_err(sck_band->sck__, 952 "%s: band %lu: invalid quadrature point index %lu.\n", 953 FUNC_NAME, (unsigned long)sck_band->id, (unsigned long)iquad_pt); 954 res = RES_BAD_ARG; 955 goto error; 956 } 957 958 quad_pt = darray_quad_pt_cdata_get(&band->quad_pts) + iquad_pt; 959 sck_quad_pt->ka_list = quad_pt->ka_list; 960 sck_quad_pt->weight = quad_pt->weight; 961 sck_quad_pt->id = iquad_pt; 962 963 exit: 964 return res; 965 error: 966 goto exit; 967 } 968 969 res_T 970 sck_band_sample_quad_pt 971 (const struct sck_band* sck_band, 972 const double r, /* Canonical random number in [0, 1[ */ 973 size_t* iquad_pt) 974 { 975 const struct band* band = NULL; 976 const double* cumul = NULL; 977 const double* find = NULL; 978 double r_next = nextafter(r, DBL_MAX); 979 size_t i; 980 981 if(!sck_band || !iquad_pt) return RES_BAD_ARG; 982 983 if(r < 0 || r >= 1) { 984 log_err(sck_band->sck__, 985 "%s: invalid canonical random number `%g'.\n", FUNC_NAME, r); 986 return RES_BAD_ARG; 987 } 988 989 band = sck_band->band__; 990 cumul = darray_double_cdata_get(&band->quad_pts_cumul); 991 992 /* Use r_next rather than r in order to find the first entry that is not less 993 * than *or equal* to r */ 994 find = search_lower_bound(&r_next, cumul, sck_band->quad_pts_count, 995 sizeof(double), cmp_dbl); 996 ASSERT(find); 997 998 /* Compute the index of the found quadrature point */ 999 i = (size_t)(find - cumul); 1000 ASSERT(i < sck_band->quad_pts_count); 1001 ASSERT(cumul[i] > r); 1002 ASSERT(!i || cumul[i-1] <= r); 1003 1004 *iquad_pt = i; 1005 1006 return RES_OK; 1007 } 1008 1009 res_T 1010 sck_quad_pt_compute_hash 1011 (const struct sck_band* band, 1012 const size_t iquad_pt, 1013 hash256_T hash) 1014 { 1015 struct sha256_ctx ctx; 1016 struct sck_quad_pt quad_pt; 1017 res_T res = RES_OK; 1018 1019 if(!band || !hash) { 1020 res = RES_BAD_ARG; 1021 goto error; 1022 } 1023 1024 res = sck_band_get_quad_pt(band, iquad_pt, &quad_pt); 1025 if(res != RES_OK) goto error; 1026 1027 sha256_ctx_init(&ctx); 1028 hash_quad_pt(&ctx, &quad_pt, band->sck__->nnodes); 1029 sha256_ctx_finalize(&ctx, hash); 1030 1031 exit: 1032 return res; 1033 error: 1034 goto exit; 1035 } 1036 1037 res_T 1038 sck_band_compute_hash 1039 (const struct sck* sck, 1040 const size_t iband, 1041 hash256_T hash) 1042 { 1043 struct sha256_ctx ctx; 1044 struct sck_band band; 1045 res_T res = RES_OK; 1046 1047 if(!sck || !hash) { 1048 res = RES_BAD_ARG; 1049 goto error; 1050 } 1051 1052 res = sck_get_band(sck, iband, &band); 1053 if(res != RES_OK) goto error; 1054 1055 sha256_ctx_init(&ctx); 1056 hash_band(&ctx, &band); 1057 sha256_ctx_finalize(&ctx, hash); 1058 1059 exit: 1060 return res; 1061 error: 1062 goto exit; 1063 } 1064 1065 res_T 1066 sck_compute_hash(const struct sck* sck, hash256_T hash) 1067 { 1068 struct sha256_ctx ctx; 1069 size_t iband; 1070 res_T res = RES_OK; 1071 1072 if(!sck || !hash) { 1073 res = RES_BAD_ARG; 1074 goto error; 1075 } 1076 1077 sha256_ctx_init(&ctx); 1078 sha256_ctx_update(&ctx, (const char*)&sck->pagesize, sizeof(sck->pagesize)); 1079 sha256_ctx_update(&ctx, (const char*)&sck->nnodes, sizeof(sck->nnodes)); 1080 FOR_EACH(iband, 0, darray_band_size_get(&sck->bands)) { 1081 struct sck_band band; 1082 SCK(get_band(sck, iband, &band)); 1083 hash_band(&ctx, &band); 1084 } 1085 sha256_ctx_finalize(&ctx, hash); 1086 1087 exit: 1088 return res; 1089 error: 1090 goto exit; 1091 } 1092 1093 const char* 1094 sck_get_name(const struct sck* sck) 1095 { 1096 ASSERT(sck); 1097 return str_cget(&sck->name); 1098 } 1099 1100 /******************************************************************************* 1101 * Local functions 1102 ******************************************************************************/ 1103 void 1104 band_release(struct band* band) 1105 { 1106 ASSERT(band); 1107 darray_quad_pt_release(&band->quad_pts); 1108 darray_double_release(&band->quad_pts_cumul); 1109 1110 if(!band->ks_list) return; 1111 1112 if(!band->map_len) { 1113 MEM_RM(band->sck->allocator, band->ks_list); 1114 } else if(band->ks_list != MAP_FAILED) { 1115 munmap(band->ks_list, band->map_len); 1116 } 1117 } 1118 1119 res_T 1120 band_copy(struct band* dst, const struct band* src) 1121 { 1122 res_T res = RES_OK; 1123 ASSERT(dst && src); 1124 1125 dst->sck = src->sck; 1126 dst->low = src->low; 1127 dst->upp = src->upp; 1128 1129 dst->map_len = src->map_len; 1130 dst->ks_list = NULL; 1131 1132 if(src->map_len) { 1133 /* The ks are mapped: copy the pointer */ 1134 dst->ks_list = src->ks_list; 1135 1136 } else if(src->ks_list != NULL) { 1137 /* The ks are loaded: duplicate table contents */ 1138 const size_t memsz = sizeof(*dst->ks_list)*src->sck->nnodes; 1139 dst->ks_list = MEM_ALLOC(dst->sck->allocator, memsz); 1140 if(!dst->ks_list) return RES_MEM_ERR; 1141 memcpy(dst->ks_list, src->ks_list, memsz); 1142 } 1143 1144 res = darray_quad_pt_copy(&dst->quad_pts, &src->quad_pts); 1145 if(res != RES_OK) return res; 1146 res = darray_double_copy(&dst->quad_pts_cumul, &src->quad_pts_cumul); 1147 if(res != RES_OK) return res; 1148 return RES_OK; 1149 } 1150 1151 void 1152 quad_pt_release(struct quad_pt* quad_pt) 1153 { 1154 ASSERT(quad_pt); 1155 if(!quad_pt->ka_list) return; 1156 1157 if(!quad_pt->map_len) { 1158 MEM_RM(quad_pt->band->sck->allocator, quad_pt->ka_list); 1159 } else if(quad_pt->ka_list != MAP_FAILED) { 1160 munmap(quad_pt->ka_list, quad_pt->map_len); 1161 } 1162 } 1163 1164 res_T 1165 quad_pt_copy(struct quad_pt* dst, const struct quad_pt* src) 1166 { 1167 ASSERT(dst && src); 1168 1169 dst->band = src->band; 1170 dst->map_len = src->map_len; 1171 dst->weight = src->weight; 1172 dst->ka_list = NULL; 1173 if(src->map_len) { 1174 /* The ka are mapped: copy the pointer */ 1175 dst->ka_list = src->ka_list; 1176 1177 } else if(src->ka_list != NULL) { 1178 /* The ka are loaded: duplicate table contents */ 1179 const size_t memsz = sizeof(*dst->ka_list)*src->band->sck->nnodes; 1180 dst->ka_list = MEM_ALLOC(dst->band->sck->allocator, memsz); 1181 if(!dst->ka_list) return RES_MEM_ERR; 1182 memcpy(dst->ka_list, src->ka_list, memsz); 1183 } 1184 return RES_OK; 1185 }