sln_stat.c (13185B)
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 file is part of Star-Line. 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #define _POSIX_C_SOURCE 200112L /* getopt */ 22 23 #include "sln.h" 24 25 #include <star/shtr.h> 26 #include <star/sbb.h> 27 #include <star/ssp.h> 28 29 #include <rsys/cstr.h> 30 #include <rsys/mem_allocator.h> 31 #include <rsys/str.h> 32 33 #include <omp.h> 34 35 #include <unistd.h> /* getopt */ 36 37 enum estimate { 38 MEAN, 39 SQMEAN, 40 ESTIMATE_COUNT__ 41 }; 42 43 #define WAVENUMBER_TO_WAVELENGTH(Nu/* [cm^-1] */) (1.e-2/(Nu))/*[m]*/ 44 45 struct args { 46 const char* tree; /* Acceleration structure */ 47 const char* molparams; 48 const char* lines; 49 50 double spectral_range[2]; /* [cm^-1]^2 */ 51 52 unsigned long nrealisations; /* Number of Monte Carlo realisations */ 53 54 /* Miscellaneous */ 55 unsigned nthreads_hint; /* Hint on the number of threads to use */ 56 int disable_line_hash_check; 57 int lines_in_shtr_format; 58 int verbose; 59 int quit; 60 }; 61 #define ARGS_DEFAULT__ {NULL,NULL,NULL,{0,DBL_MAX},10000,UINT_MAX,0,0,0,0} 62 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 63 64 struct cmd { 65 struct args args; 66 67 struct sln_tree* tree; 68 unsigned nthreads; 69 }; 70 #define CMD_NULL__ {0} 71 static const struct cmd CMD_NULL = CMD_NULL__; 72 73 struct accum { 74 double sum; 75 double sum2; 76 size_t count; 77 }; 78 #define ACCUM_NULL__ {0} 79 80 /******************************************************************************* 81 * Helper functions 82 ******************************************************************************/ 83 static void 84 usage(FILE* stream) 85 { 86 fprintf(stream, 87 "usage: sln-stat [-dhsv] [-n nrealisations] [-t threads]\n" 88 " -S nu_min,nu_max -a accel_struct -m molparams -l lines\n"); 89 } 90 91 static res_T 92 parse_spectral_range(const char* str, double spectral_range[2]) 93 { 94 size_t len = 0; 95 res_T res = RES_OK; 96 ASSERT(str && spectral_range); 97 98 res = cstr_to_list_double(str, ',', spectral_range, &len, 2); 99 if(res == RES_OK && len < 2) res = RES_BAD_ARG; 100 101 return res; 102 } 103 104 static res_T 105 args_init(struct args* args, int argc, char** argv) 106 { 107 int opt = 0; 108 res_T res = RES_OK; 109 110 ASSERT(args); 111 112 *args = ARGS_DEFAULT; 113 114 while((opt = getopt(argc, argv, "a:dhl:m:n:S:st:v")) != -1) { 115 switch(opt) { 116 case 'a': args->tree = optarg; break; 117 case 'd': args->disable_line_hash_check = 1; break; 118 case 'h': 119 usage(stdout); 120 args->quit = 1; 121 goto exit; 122 case 'l': args->lines = optarg; break; 123 case 'm': args->molparams = optarg; break; 124 case 'n': res = cstr_to_ulong(optarg, &args->nrealisations); break; 125 case 'S': res = parse_spectral_range(optarg, args->spectral_range); break; 126 case 's': args->lines_in_shtr_format = 1; break; 127 case 't': 128 res = cstr_to_uint(optarg, &args->nthreads_hint); 129 if(res == RES_OK && args->nthreads_hint == 0) res = RES_BAD_ARG; 130 break; 131 case 'v': args->verbose += (args->verbose < 3); break; 132 default: res = RES_BAD_ARG; break; 133 } 134 if(res != RES_OK) { 135 if(optarg) { 136 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 137 argv[0], optarg, opt); 138 } 139 goto error; 140 } 141 } 142 143 #define MANDATORY(Cond, Name, Opt) { \ 144 if(!(Cond)) { \ 145 fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \ 146 res = RES_BAD_ARG; \ 147 goto error; \ 148 } \ 149 } (void)0 150 MANDATORY(args->molparams, "molparams", 'm'); 151 MANDATORY(args->lines, "line list", 'l'); 152 MANDATORY(args->tree, "acceleration structure", 'a'); 153 #undef MANDATORY 154 155 exit: 156 return res; 157 error: 158 usage(stderr); 159 goto exit; 160 } 161 162 static res_T 163 load_lines 164 (struct shtr* shtr, 165 const struct args* args, 166 struct shtr_line_list** out_lines) 167 { 168 struct shtr_line_list* lines = NULL; 169 res_T res = RES_OK; 170 ASSERT(shtr && args && out_lines); 171 172 if(args->lines_in_shtr_format) { 173 struct shtr_line_list_read_args read_args = SHTR_LINE_LIST_READ_ARGS_NULL; 174 175 /* Loads lines from data serialized by the Star-HITRAN library */ 176 read_args.filename = args->lines; 177 res = shtr_line_list_read(shtr, &read_args, &lines); 178 if(res != RES_OK) goto error; 179 180 } else { 181 struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL; 182 183 /* Loads lines from a file in HITRAN format */ 184 load_args.filename = args->lines; 185 res = shtr_line_list_load(shtr, &load_args, &lines); 186 if(res != RES_OK) goto error; 187 } 188 189 exit: 190 *out_lines = lines; 191 return res; 192 error: 193 if(lines) { SHTR(line_list_ref_put(lines)); lines = NULL; } 194 goto exit; 195 } 196 197 static void 198 delete_per_thread_rngs(const struct cmd* cmd, struct ssp_rng* rngs[]) 199 { 200 unsigned i = 0; 201 ASSERT(cmd && rngs); 202 203 FOR_EACH(i, 0, cmd->nthreads) { 204 if(rngs[i]) SSP(rng_ref_put(rngs[i])); 205 } 206 mem_rm(rngs); 207 } 208 209 static res_T 210 create_per_thread_rngs(const struct cmd* cmd, struct ssp_rng** out_rngs[]) 211 { 212 struct ssp_rng_proxy* proxy = NULL; 213 struct ssp_rng** rngs = NULL; 214 size_t i = 0; 215 res_T res = RES_OK; 216 ASSERT(cmd); 217 218 rngs = mem_calloc(cmd->nthreads, sizeof(*rngs)); 219 if(!rngs) { res = RES_MEM_ERR; goto error; } 220 221 res = ssp_rng_proxy_create(NULL, SSP_RNG_THREEFRY, cmd->nthreads, &proxy); 222 if(res != RES_OK) goto error; 223 224 FOR_EACH(i, 0, cmd->nthreads) { 225 res = ssp_rng_proxy_create_rng(proxy, i, &rngs[i]); 226 if(res != RES_OK) goto error; 227 } 228 229 exit: 230 *out_rngs = rngs; 231 if(proxy) SSP(rng_proxy_ref_put(proxy)); 232 return res; 233 error: 234 if(cmd->args.verbose >= 1) { 235 fprintf(stderr, 236 "Error creating the list of per thread RNG -- %s\n", 237 res_to_cstr(res)); 238 } 239 if(rngs) delete_per_thread_rngs(cmd, rngs); 240 rngs = NULL; 241 goto exit; 242 } 243 244 static void 245 cmd_release(struct cmd* cmd) 246 { 247 ASSERT(cmd); 248 if(cmd->tree) SLN(tree_ref_put(cmd->tree)); 249 } 250 251 static res_T 252 cmd_init(struct cmd* cmd, const struct args* args) 253 { 254 /* Star Line */ 255 struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 256 struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL; 257 struct sln_device* sln = NULL; 258 259 /* Star HITRAN */ 260 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 261 struct shtr* shtr = NULL; 262 struct shtr_isotope_metadata* molparams = NULL; 263 struct shtr_line_list* lines = NULL; 264 265 /* Miscellaneous */ 266 unsigned nthreads_max = 0; 267 res_T res = RES_OK; 268 269 ASSERT(cmd && args); 270 271 *cmd = CMD_NULL; 272 273 shtr_args.verbose = args->verbose; 274 res = shtr_create(&shtr_args, &shtr); 275 if(res != RES_OK) goto error; 276 277 res = shtr_isotope_metadata_load(shtr, args->molparams, &molparams); 278 if(res != RES_OK) goto error; 279 280 res = load_lines(shtr, args, &lines); 281 if(res != RES_OK) goto error; 282 283 sln_args.verbose = args->verbose; 284 res = sln_device_create(&sln_args, &sln); 285 if(res != RES_OK) goto error; 286 287 tree_args.metadata = molparams; 288 tree_args.lines = lines; 289 tree_args.filename = args->tree; 290 tree_args.disable_line_hash_check = args->disable_line_hash_check; 291 res = sln_tree_read(sln, &tree_args, &cmd->tree); 292 if(res != RES_OK) goto error; 293 294 nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs()); 295 cmd->args = *args; 296 cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max); 297 298 exit: 299 if(sln) SLN(device_ref_put(sln)); 300 if(shtr) SHTR(ref_put(shtr)); 301 if(molparams) SHTR(isotope_metadata_ref_put(molparams)); 302 if(lines) SHTR(line_list_ref_put(lines)); 303 return res; 304 error: 305 cmd_release(cmd); 306 *cmd = CMD_NULL; 307 goto exit; 308 } 309 310 static INLINE const char* 311 estimate_cstr(const enum estimate estimate) 312 { 313 const char* cstr = NULL; 314 switch(estimate) { 315 case MEAN: cstr="mean"; break; 316 case SQMEAN: cstr="mean_of_squares"; break; 317 default: FATAL("Unreachable code\n"); break; 318 } 319 return cstr; 320 } 321 322 static res_T 323 realisation 324 (const struct cmd* cmd, 325 struct ssp_rng* rng, 326 double out_weights[ESTIMATE_COUNT__]) 327 { 328 /* Acceleration structure */ 329 struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL; 330 const struct sln_node* root = NULL; 331 332 /* Variables to sample k */ 333 const struct sln_node* leaf1 = NULL; 334 const struct sln_node* leaf2 = NULL; 335 double leaf_proba1 = 0; /* Probability of sampling a line */ 336 double leaf_proba2 = 0; /* Probability of sampling a line */ 337 double leaf_ka1 = 0; /* Value of a line */ 338 double leaf_ka2 = 0; /* Value of a line */ 339 340 /* Miscellaneous */ 341 double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */ 342 double nu = 0; /* Sampled wavenumber [cm^-1] */ 343 int i = 0; 344 res_T res = RES_OK; 345 346 ASSERT(cmd && rng && out_weights); /* Check pre-conditions */ 347 348 /* Uniformly sample the spectral dimension */ 349 nu = ssp_rng_uniform_double 350 (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]); 351 352 SLN(tree_get_desc(cmd->tree, &tree_desc)); 353 354 /* Store the root node of the tree */ 355 root = sln_tree_get_root(cmd->tree); 356 357 /* Importance sampling of a line and evaluation of the line contribution */ 358 leaf1 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba1); 359 if(!leaf1) { res = RES_BAD_ARG; goto error; } 360 leaf_ka1 = sln_node_eval(cmd->tree, leaf1, NULL, nu); 361 362 /* Importance sampling of a line and evaluation of the line contribution */ 363 leaf2 = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba2); 364 if(!leaf2) { res = RES_BAD_ARG; goto error; } 365 leaf_ka2 = sln_node_eval(cmd->tree, leaf2, NULL, nu); 366 367 w[MEAN] = leaf_ka1 / leaf_proba1; 368 w[SQMEAN] = (leaf_ka1/leaf_proba1)*(leaf_ka2/leaf_proba2); 369 370 exit: 371 FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i]; 372 return res; 373 error: 374 FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN; 375 goto exit; 376 } 377 378 static res_T 379 cmd_run(const struct cmd* cmd) 380 { 381 /* Random Number Generator */ 382 struct ssp_rng** rngs = NULL; 383 384 /* Monte Carlo */ 385 struct accum accum[ESTIMATE_COUNT__] = {0}; 386 int64_t i = 0; /* Index of the realisation */ 387 size_t nrejects = 0; /* Number of rejected realisations */ 388 389 /* Progress */ 390 size_t nrealisations = 0; 391 size_t realisation_done = 0; 392 int progress = 0; 393 int progress_pcent = 10; 394 395 res_T res = RES_OK; 396 ASSERT(cmd); 397 398 res = create_per_thread_rngs(cmd, &rngs); 399 if(res != RES_OK) goto error; 400 401 #define PROGRESS_MSG "Solving: %3d%%\n" 402 if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress); 403 404 nrealisations = cmd->args.nrealisations; 405 406 omp_set_num_threads((int)cmd->nthreads); 407 408 #pragma omp parallel for schedule(static) 409 for(i = 0; i < (int64_t)nrealisations; ++i) { 410 double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */ 411 const int ithread = omp_get_thread_num(); 412 int pcent = 0; 413 res_T res_realisation = RES_OK; 414 415 res_realisation = realisation(cmd, rngs[ithread], w); 416 417 #pragma omp critical 418 { 419 /* Update the Monte Carlo accumulator */ 420 if(res_realisation == RES_OK) { 421 int iestim = 0; 422 FOR_EACH(iestim, 0, ESTIMATE_COUNT__) { 423 accum[iestim].sum += w[iestim]; 424 accum[iestim].sum2 += w[iestim]*w[iestim]; 425 accum[iestim].count += 1; 426 } 427 } 428 429 if(cmd->args.verbose >= 3) { 430 /* Update progress */ 431 realisation_done += 1; 432 pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5); 433 if(pcent/progress_pcent > progress/progress_pcent) { 434 progress = pcent; 435 fprintf(stderr, PROGRESS_MSG, progress); 436 } 437 } 438 } 439 } 440 441 #undef PROGRESS_MSG 442 443 nrejects = nrealisations - accum[0].count; 444 445 FOR_EACH(i, 0, ESTIMATE_COUNT__) { 446 const double E = accum[i].sum / (double)accum[i].count; 447 const double V = accum[i].sum2 / (double)accum[i].count - E*E; 448 const double SE = sqrt(V/(double)accum[i].count); 449 450 /* Assume that the number of realisations is the same for all estimates */ 451 ASSERT(accum[i].count == accum[0].count); 452 453 printf("%-16s: %e +/- %e; %lu\n", 454 estimate_cstr(i), E, SE, (unsigned long)nrejects); 455 } 456 457 exit: 458 delete_per_thread_rngs(cmd, rngs); 459 return res; 460 error: 461 goto exit; 462 } 463 464 /******************************************************************************* 465 * Main function 466 ******************************************************************************/ 467 int 468 main(int argc, char** argv) 469 { 470 struct args args = ARGS_DEFAULT; 471 struct cmd cmd = CMD_NULL; 472 int err = 0; 473 res_T res = RES_OK; 474 475 if((res = args_init(&args, argc, argv)) != RES_OK) goto error; 476 if(args.quit) goto exit; 477 478 if((res = cmd_init(&cmd, &args)) != RES_OK) goto error; 479 if((res = cmd_run(&cmd)) != RES_OK) goto error; 480 481 exit: 482 cmd_release(&cmd); 483 CHK(mem_allocated_size() == 0); 484 return err; 485 error: 486 err = 1; 487 goto exit; 488 }