ssol_estimator_c.h (19284B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 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 #ifndef SSOL_ESTIMATOR_C_H 18 #define SSOL_ESTIMATOR_C_H 19 20 #include "ssol_device_c.h" 21 #include "ssol_instance_c.h" 22 #include "ssol_shape_c.h" 23 24 #include <rsys/ref_count.h> 25 #include <rsys/hash_table.h> 26 27 /* Forward declaration */ 28 struct mem_allocator; 29 struct ssol_instance; 30 struct ssp_rng_proxy; 31 32 /* Monte carlo data */ 33 struct mc_data { 34 size_t irealisation; 35 double tmp; 36 37 /* Internal data; use get() */ 38 double weight__; 39 double sqr_weight__; 40 }; 41 #define MC_DATA_NULL__ { SIZE_MAX, 0, 0, 0 } 42 static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__; 43 44 #define MC_RECEIVER_DATA \ 45 struct mc_data incoming_flux; /* In W */ \ 46 struct mc_data incoming_if_no_atm_loss; /* In W */ \ 47 struct mc_data incoming_if_no_field_loss; /* In W */ \ 48 struct mc_data incoming_lost_in_field; /* In W */ \ 49 struct mc_data incoming_lost_in_atmosphere; /* In W */ \ 50 struct mc_data absorbed_flux; /* In W */ \ 51 struct mc_data absorbed_if_no_atm_loss; /* In W */ \ 52 struct mc_data absorbed_if_no_field_loss; /* In W */ \ 53 struct mc_data absorbed_lost_in_field; /* In W */ \ 54 struct mc_data absorbed_lost_in_atmosphere; /* In W */ 55 56 /******************************************************************************* 57 * Deferred MC data accumulators 58 ******************************************************************************/ 59 static INLINE void 60 mc_data_init(struct mc_data* data) 61 { 62 ASSERT(data); 63 *data = MC_DATA_NULL; 64 } 65 66 static INLINE void 67 mc_data_flush(struct mc_data* data) 68 { 69 ASSERT(data); 70 data->weight__ += data->tmp; 71 data->sqr_weight__ += data->tmp * data->tmp; 72 data->tmp = 0; 73 } 74 75 static INLINE void 76 mc_data_add_weight 77 (struct mc_data* data, 78 const size_t irealisation, 79 const double w) 80 { 81 ASSERT(data); 82 ASSERT(irealisation != SIZE_MAX); 83 if(irealisation != data->irealisation) { 84 mc_data_flush(data); 85 data->irealisation = irealisation; 86 } 87 data->tmp += w; 88 } 89 90 static INLINE void 91 mc_data_accum(struct mc_data* dst, struct mc_data* src) 92 { 93 ASSERT(dst && src); 94 mc_data_flush(dst); 95 mc_data_flush(src); 96 dst->weight__ += src->weight__; 97 dst->sqr_weight__ += src->sqr_weight__; 98 } 99 100 static INLINE void 101 mc_data_get(struct mc_data* data, double* weight, double* sqr_weight) 102 { 103 ASSERT(data && weight && sqr_weight); 104 mc_data_flush(data); 105 *weight = data->weight__; 106 *sqr_weight = data->sqr_weight__; 107 } 108 109 static INLINE void 110 mc_data_cancel(struct mc_data* data, const size_t irealisation) 111 { 112 ASSERT(data); 113 if(data->irealisation!=irealisation) return; 114 data->tmp = 0; 115 } 116 117 static INLINE void 118 mc_data_apply_factor 119 (struct mc_data* data, 120 const size_t irealisation, 121 const double factor) 122 { 123 ASSERT(data); 124 if(data->irealisation != irealisation) return; 125 data->tmp *= factor; 126 } 127 128 /******************************************************************************* 129 * One sided per shape MC data 130 ******************************************************************************/ 131 struct mc_primitive_1side { 132 MC_RECEIVER_DATA 133 }; 134 135 #define MC_PRIMITIVE_1SIDE_NULL__ { \ 136 MC_DATA_NULL__, \ 137 MC_DATA_NULL__, \ 138 MC_DATA_NULL__, \ 139 MC_DATA_NULL__, \ 140 MC_DATA_NULL__, \ 141 MC_DATA_NULL__, \ 142 MC_DATA_NULL__, \ 143 MC_DATA_NULL__, \ 144 MC_DATA_NULL__, \ 145 MC_DATA_NULL__ \ 146 } 147 148 static const struct mc_primitive_1side MC_PRIMITIVE_1SIDE_NULL = 149 MC_PRIMITIVE_1SIDE_NULL__; 150 151 /* Map an unsigned to a struct mc_primitive_1side */ 152 #define HTABLE_NAME prim2mc 153 #define HTABLE_KEY unsigned 154 #define HTABLE_DATA struct mc_primitive_1side 155 #include <rsys/hash_table.h> 156 157 struct mc_shape_1side { 158 struct htable_prim2mc prim2mc; 159 }; 160 161 static INLINE void 162 mc_shape_1side_init 163 (struct mem_allocator* allocator, struct mc_shape_1side* mc) 164 { 165 ASSERT(mc); 166 htable_prim2mc_init(allocator, &mc->prim2mc); 167 } 168 169 static INLINE void 170 mc_shape_1side_release(struct mc_shape_1side* mc) 171 { 172 ASSERT(mc); 173 htable_prim2mc_release(&mc->prim2mc); 174 } 175 176 static INLINE res_T 177 mc_shape_1side_copy 178 (struct mc_shape_1side* dst, const struct mc_shape_1side* src) 179 { 180 ASSERT(dst && src); 181 return htable_prim2mc_copy(&dst->prim2mc, &src->prim2mc); 182 } 183 184 static INLINE res_T 185 mc_shape_1side_copy_and_release 186 (struct mc_shape_1side* dst, struct mc_shape_1side* src) 187 { 188 ASSERT(dst && src); 189 return htable_prim2mc_copy_and_release(&dst->prim2mc, &src->prim2mc); 190 } 191 192 static INLINE res_T 193 mc_shape_1side_get_mc_primitive 194 (struct mc_shape_1side* mc_shape1, 195 const unsigned iprim, 196 struct mc_primitive_1side** out_mc_prim1) 197 { 198 struct mc_primitive_1side* mc_prim1 = NULL; 199 res_T res = RES_OK; 200 ASSERT(mc_shape1 && out_mc_prim1); 201 202 mc_prim1 = htable_prim2mc_find(&mc_shape1->prim2mc, &iprim); 203 if(!mc_prim1) { 204 res = htable_prim2mc_set(&mc_shape1->prim2mc, &iprim, &MC_PRIMITIVE_1SIDE_NULL); 205 if(res != RES_OK) goto error; 206 207 mc_prim1 = htable_prim2mc_find(&mc_shape1->prim2mc, &iprim); 208 } 209 210 exit: 211 *out_mc_prim1 = mc_prim1; 212 return res; 213 error: 214 goto exit; 215 } 216 217 /******************************************************************************* 218 * One sided per receiver MC data 219 ******************************************************************************/ 220 /* Map a ssol shape to a struct mc_shape_1side */ 221 #define HTABLE_NAME shape2mc 222 #define HTABLE_KEY const struct ssol_shape* 223 #define HTABLE_DATA struct mc_shape_1side 224 #define HTABLE_DATA_FUNCTOR_INIT mc_shape_1side_init 225 #define HTABLE_DATA_FUNCTOR_RELEASE mc_shape_1side_release 226 #define HTABLE_DATA_FUNCTOR_COPY mc_shape_1side_copy 227 #define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE mc_shape_1side_copy_and_release 228 #include <rsys/hash_table.h> 229 230 struct mc_receiver_1side { 231 MC_RECEIVER_DATA 232 struct htable_shape2mc shape2mc; 233 }; 234 235 static FINLINE void 236 mc_receiver_1side_copy_mc_weights__ 237 (struct mc_receiver_1side* dst, const struct mc_receiver_1side* src) 238 { 239 ASSERT(dst && src); 240 dst->incoming_flux = src->incoming_flux; 241 dst->incoming_if_no_atm_loss = src->incoming_if_no_atm_loss; 242 dst->incoming_if_no_field_loss = src->incoming_if_no_field_loss; 243 dst->incoming_lost_in_atmosphere = src->incoming_lost_in_atmosphere; 244 dst->incoming_lost_in_field = src->incoming_lost_in_field; 245 dst->absorbed_flux = src->absorbed_flux; 246 dst->absorbed_if_no_atm_loss = src->absorbed_if_no_atm_loss; 247 dst->absorbed_if_no_field_loss = src->absorbed_if_no_field_loss; 248 dst->absorbed_lost_in_atmosphere = src->absorbed_lost_in_atmosphere; 249 dst->absorbed_lost_in_field = src->absorbed_lost_in_field; 250 } 251 252 static INLINE void 253 mc_receiver_1side_init 254 (struct mem_allocator* allocator, struct mc_receiver_1side* mc) 255 { 256 ASSERT(mc); 257 mc->incoming_flux = MC_DATA_NULL; 258 mc->incoming_if_no_atm_loss = MC_DATA_NULL; 259 mc->incoming_if_no_field_loss = MC_DATA_NULL; 260 mc->incoming_lost_in_atmosphere = MC_DATA_NULL; 261 mc->incoming_lost_in_field = MC_DATA_NULL; 262 mc->absorbed_flux = MC_DATA_NULL; 263 mc->absorbed_if_no_atm_loss = MC_DATA_NULL; 264 mc->absorbed_if_no_field_loss = MC_DATA_NULL; 265 mc->absorbed_lost_in_atmosphere = MC_DATA_NULL; 266 mc->absorbed_lost_in_field = MC_DATA_NULL; 267 htable_shape2mc_init(allocator, &mc->shape2mc); 268 } 269 270 static INLINE void 271 mc_receiver_1side_release(struct mc_receiver_1side* mc) 272 { 273 ASSERT(mc); 274 htable_shape2mc_release(&mc->shape2mc); 275 } 276 277 static INLINE res_T 278 mc_receiver_1side_copy 279 (struct mc_receiver_1side* dst, const struct mc_receiver_1side* src) 280 { 281 ASSERT(dst && src); 282 mc_receiver_1side_copy_mc_weights__(dst, src); 283 return htable_shape2mc_copy(&dst->shape2mc, &src->shape2mc); 284 } 285 286 static INLINE res_T 287 mc_receiver_1side_copy_and_release 288 (struct mc_receiver_1side* dst, struct mc_receiver_1side* src) 289 { 290 ASSERT(dst && src); 291 mc_receiver_1side_copy_mc_weights__(dst, src); 292 return htable_shape2mc_copy_and_release(&dst->shape2mc, &src->shape2mc); 293 } 294 295 static INLINE res_T 296 mc_receiver_1side_get_mc_shape 297 (struct mc_receiver_1side* mc_rcv, 298 const struct ssol_shape* shape, 299 struct mc_shape_1side** out_mc_shape1) 300 { 301 struct mc_shape_1side* mc_shape1 = NULL; 302 struct mc_shape_1side mc_shape1_null; 303 res_T res = RES_OK; 304 ASSERT(mc_rcv && shape && out_mc_shape1); 305 306 mc_shape_1side_init(shape->dev->allocator, &mc_shape1_null); 307 308 mc_shape1 = htable_shape2mc_find(&mc_rcv->shape2mc, &shape); 309 if(!mc_shape1) { 310 res = htable_shape2mc_set(&mc_rcv->shape2mc, &shape, &mc_shape1_null); 311 if(res != RES_OK) goto error; 312 313 mc_shape1 = htable_shape2mc_find(&mc_rcv->shape2mc, &shape); 314 } 315 316 exit: 317 mc_shape_1side_release(&mc_shape1_null); 318 *out_mc_shape1 = mc_shape1; 319 return res; 320 error: 321 goto exit; 322 } 323 324 /******************************************************************************* 325 * Double sided per receiver MC data 326 ******************************************************************************/ 327 struct mc_receiver { 328 struct mc_receiver_1side front; 329 struct mc_receiver_1side back; 330 }; 331 332 static INLINE void 333 mc_receiver_init(struct mem_allocator* allocator, struct mc_receiver* mc) 334 { 335 ASSERT(mc); 336 mc_receiver_1side_init(allocator, &mc->front); 337 mc_receiver_1side_init(allocator, &mc->back); 338 } 339 340 static INLINE void 341 mc_receiver_release(struct mc_receiver* mc) 342 { 343 ASSERT(mc); 344 mc_receiver_1side_release(&mc->front); 345 mc_receiver_1side_release(&mc->back); 346 } 347 348 static INLINE res_T 349 mc_receiver_copy(struct mc_receiver* dst, const struct mc_receiver* src) 350 { 351 res_T res = RES_OK; 352 ASSERT(dst && src); 353 res = mc_receiver_1side_copy(&dst->front, &src->front); 354 if(res != RES_OK) return res; 355 res = mc_receiver_1side_copy(&dst->back, &src->back); 356 if(res != RES_OK) return res; 357 return RES_OK; 358 } 359 360 static INLINE res_T 361 mc_receiver_copy_and_release 362 (struct mc_receiver* dst, struct mc_receiver* src) 363 { 364 res_T res = RES_OK; 365 ASSERT(dst && src); 366 res = mc_receiver_1side_copy_and_release(&dst->front, &src->front); 367 if(res != RES_OK) return res; 368 res = mc_receiver_1side_copy_and_release(&dst->back, &src->back); 369 if(res != RES_OK) return res; 370 return RES_OK; 371 } 372 373 /* Define the htable_receiver data structure */ 374 #define HTABLE_NAME receiver 375 #define HTABLE_KEY const struct ssol_instance* 376 #define HTABLE_DATA struct mc_receiver 377 #define HTABLE_DATA_FUNCTOR_INIT mc_receiver_init 378 #define HTABLE_DATA_FUNCTOR_RELEASE mc_receiver_release 379 #define HTABLE_DATA_FUNCTOR_COPY mc_receiver_copy 380 #define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE mc_receiver_copy_and_release 381 #include <rsys/hash_table.h> 382 383 /******************************************************************************* 384 * Per sampled instance MC data 385 ******************************************************************************/ 386 struct mc_sampled { 387 /* Global data for this entity */ 388 struct mc_data cos_factor; 389 struct mc_data shadowed; 390 size_t nb_samples; 391 392 /* By-receptor data for this entity */ 393 struct htable_receiver mc_rcvs; 394 }; 395 396 static INLINE void 397 mc_sampled_init 398 (struct mem_allocator* allocator, 399 struct mc_sampled* samp) 400 { 401 ASSERT(samp); 402 samp->cos_factor = MC_DATA_NULL; 403 samp->shadowed = MC_DATA_NULL; 404 samp->nb_samples = 0; 405 htable_receiver_init(allocator, &samp->mc_rcvs); 406 } 407 408 static INLINE void 409 mc_sampled_release(struct mc_sampled* samp) 410 { 411 ASSERT(samp); 412 htable_receiver_release(&samp->mc_rcvs); 413 } 414 415 static INLINE res_T 416 mc_sampled_copy(struct mc_sampled* dst, const struct mc_sampled* src) 417 { 418 ASSERT(dst && src); 419 dst->cos_factor = src->cos_factor; 420 dst->shadowed = src->shadowed; 421 dst->nb_samples = src->nb_samples; 422 return htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs); 423 } 424 425 static INLINE res_T 426 mc_sampled_copy_and_release(struct mc_sampled* dst, struct mc_sampled* src) 427 { 428 ASSERT(dst && src); 429 dst->cos_factor = src->cos_factor; 430 dst->shadowed = src->shadowed; 431 dst->nb_samples = src->nb_samples; 432 return htable_receiver_copy_and_release(&dst->mc_rcvs, &src->mc_rcvs); 433 } 434 435 static INLINE res_T 436 mc_sampled_get_mc_receiver_1side 437 (struct mc_sampled* mc_samp, 438 const struct ssol_instance* inst, 439 const enum ssol_side_flag side, 440 struct mc_receiver_1side** out_mc_rcv1) 441 { 442 struct mc_receiver* mc_rcv = NULL; 443 struct mc_receiver_1side* mc_rcv1 = NULL; 444 struct mc_receiver mc_rcv_null; 445 res_T res = RES_OK; 446 ASSERT(mc_samp && inst); 447 ASSERT(inst->receiver_mask & (int)side); 448 449 mc_receiver_init(inst->dev->allocator, &mc_rcv_null); 450 451 mc_rcv = htable_receiver_find(&mc_samp->mc_rcvs, &inst); 452 if(!mc_rcv) { 453 res = htable_receiver_set(&mc_samp->mc_rcvs, &inst, &mc_rcv_null); 454 if(res != RES_OK) goto error; 455 mc_rcv = htable_receiver_find(&mc_samp->mc_rcvs, &inst); 456 } 457 mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; 458 459 exit: 460 mc_receiver_release(&mc_rcv_null); 461 *out_mc_rcv1 = mc_rcv1; 462 return res; 463 error: 464 mc_rcv1 = NULL; 465 goto exit; 466 } 467 468 /* Define the htable_primary data structure */ 469 #define HTABLE_NAME sampled 470 #define HTABLE_KEY const struct ssol_instance* 471 #define HTABLE_DATA struct mc_sampled 472 #define HTABLE_DATA_FUNCTOR_INIT mc_sampled_init 473 #define HTABLE_DATA_FUNCTOR_RELEASE mc_sampled_release 474 #define HTABLE_DATA_FUNCTOR_COPY mc_sampled_copy 475 #define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE mc_sampled_copy_and_release 476 #include <rsys/hash_table.h> 477 478 /******************************************************************************* 479 * Radiative path 480 ******************************************************************************/ 481 #define DARRAY_NAME path_vertex 482 #define DARRAY_DATA struct ssol_path_vertex 483 #include <rsys/dynamic_array.h> 484 485 struct path { 486 enum ssol_path_type type; 487 struct darray_path_vertex vertices; 488 }; 489 490 static INLINE void 491 path_init(struct mem_allocator* allocator, struct path* path) 492 { 493 ASSERT(path); 494 path->type = SSOL_PATH_MISSING; 495 darray_path_vertex_init(allocator, &path->vertices); 496 } 497 498 static INLINE void 499 path_release(struct path* path) 500 { 501 ASSERT(path); 502 darray_path_vertex_release(&path->vertices); 503 } 504 505 static INLINE res_T 506 path_copy(struct path* dst, const struct path* src) 507 { 508 ASSERT(dst && src); 509 dst->type = src->type; 510 return darray_path_vertex_copy(&dst->vertices, &src->vertices); 511 } 512 513 static INLINE res_T 514 path_copy_and_release(struct path* dst, struct path* src) 515 { 516 ASSERT(dst && src); 517 dst->type = src->type; 518 return darray_path_vertex_copy_and_release(&dst->vertices, &src->vertices); 519 } 520 521 static INLINE res_T 522 path_copy_and_clear(struct path* dst, struct path* src) 523 { 524 ASSERT(dst && src); 525 dst->type = src->type; 526 return darray_path_vertex_copy_and_clear(&dst->vertices, &src->vertices); 527 } 528 529 static INLINE res_T 530 path_add_vertex(struct path* path, const double pos[3], const double weight) 531 { 532 struct ssol_path_vertex vertex; 533 ASSERT(path && pos && weight >= 0); 534 vertex.pos[0] = pos[0]; 535 vertex.pos[1] = pos[1]; 536 vertex.pos[2] = pos[2]; 537 vertex.weight = weight; 538 return darray_path_vertex_push_back(&path->vertices, &vertex); 539 } 540 541 #define DARRAY_NAME path 542 #define DARRAY_DATA struct path 543 #define DARRAY_FUNCTOR_INIT path_init 544 #define DARRAY_FUNCTOR_RELEASE path_release 545 #define DARRAY_FUNCTOR_COPY path_copy 546 #define DARRAY_FUNCTOR_COPY_AND_RELEASE path_copy_and_release 547 #include <rsys/dynamic_array.h> 548 549 /******************************************************************************* 550 * Estimator data structure 551 ******************************************************************************/ 552 struct ssol_estimator { 553 size_t realisation_count; 554 size_t failed_count; 555 556 /* Implicit MC computations */ 557 struct mc_data cos_factor; 558 struct mc_data absorbed_by_receivers; 559 struct mc_data shadowed; 560 struct mc_data missing; 561 struct mc_data extinguished_by_atmosphere; 562 struct mc_data other_absorbed; 563 564 struct htable_receiver mc_receivers; /* Per receiver MC */ 565 struct htable_sampled mc_sampled; /* Per sampled instance MC */ 566 567 struct darray_path paths; /* Tracked paths */ 568 569 /* Overall area of the sampled instances. Actually this is not the area that 570 * is effectively sampled since an instance may be sampled through a proxy 571 * geometry */ 572 double sampled_area; 573 574 /* State of the RNG after the simulation */ 575 struct ssp_rng* rng; 576 577 struct ssol_device* dev; 578 ref_T ref; 579 }; 580 581 extern LOCAL_SYM res_T 582 estimator_create 583 (struct ssol_device* dev, 584 struct ssol_scene* scene, 585 struct ssol_estimator** estimator); 586 587 extern LOCAL_SYM res_T 588 estimator_save_rng_state 589 (struct ssol_estimator* estimator, 590 const struct ssp_rng_proxy* proxy); 591 592 static FINLINE res_T 593 get_mc_receiver_1side 594 (struct htable_receiver* receivers, 595 const struct ssol_instance* inst, 596 const enum ssol_side_flag side, 597 struct mc_receiver_1side** out_mc_rcv1) 598 { 599 struct mc_receiver* mc_rcv = NULL; 600 struct mc_receiver_1side* mc_rcv1 = NULL; 601 struct mc_receiver mc_rcv_null; 602 res_T res = RES_OK; 603 ASSERT(receivers && inst && out_mc_rcv1); 604 ASSERT(inst->receiver_mask & (int)side); 605 606 mc_receiver_init(inst->dev->allocator, &mc_rcv_null); 607 608 mc_rcv = htable_receiver_find(receivers, &inst); 609 if(!mc_rcv) { 610 res = htable_receiver_set(receivers, &inst, &mc_rcv_null); 611 if(res != RES_OK) goto error; 612 mc_rcv = htable_receiver_find(receivers, &inst); 613 } 614 615 mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back; 616 exit: 617 mc_receiver_release(&mc_rcv_null); 618 *out_mc_rcv1 = mc_rcv1; 619 return res; 620 error: 621 goto exit; 622 } 623 624 static FINLINE res_T 625 get_mc_sampled 626 (struct htable_sampled* sampled, 627 const struct ssol_instance* inst, 628 struct mc_sampled** out_mc_samp) 629 { 630 struct mc_sampled* mc_samp = NULL; 631 res_T res = RES_OK; 632 ASSERT(sampled && inst && out_mc_samp); 633 634 mc_samp = htable_sampled_find(sampled, &inst); 635 if(!mc_samp) { 636 struct mc_sampled mc_samp_null; 637 mc_sampled_init(inst->dev->allocator, &mc_samp_null); 638 res = htable_sampled_set(sampled, &inst, &mc_samp_null); 639 mc_sampled_release(&mc_samp_null); 640 if(res != RES_OK) goto error; 641 mc_samp = htable_sampled_find(sampled, &inst); 642 } 643 644 exit: 645 *out_mc_samp = mc_samp; 646 return res; 647 error: 648 goto exit; 649 } 650 651 #endif /* SSOL_ESTIMATOR_C_H */ 652