test_solstice_simulation.c (14877B)
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 #define _POSIX_C_SOURCE 200809L /* mkstemp support */ 18 19 #include <rsys/rsys.h> 20 #include <rsys/math.h> 21 #include <rsys/double2.h> 22 23 #include <math.h> 24 #include <stdio.h> 25 #include <stdlib.h> 26 #include <string.h> 27 28 #ifdef COMPILER_CL 29 /* Wrap POSIX functions and constants */ 30 #include <io.h> 31 #include <fcntl.h> 32 #include <sys/stat.h> 33 #define fdopen _fdopen 34 #define open _open 35 #endif 36 37 enum side { 38 FRONT, 39 BACK 40 }; 41 42 enum global_result_type { 43 GLOBAL_POTENTIAL, 44 GLOBAL_ABSORBED, 45 GLOBAL_COS, 46 GLOBAL_SHADOW, 47 GLOBAL_MISSING, 48 GLOBAL_ATMOSPHERE, 49 GLOBAL_REFLECTIVITY, 50 GLOBAL_RESULTS_COUNT__ 51 }; 52 53 enum receiver_result_type { 54 FIRST_RECEIVER_RESULT, 55 FRONT_ABSORBED_FLUX = FIRST_RECEIVER_RESULT, 56 FRONT_INCOMING_FLUX, 57 FRONT_ABSORBED_FIELD_GAIN, 58 FRONT_ABSORBED_ATM_GAIN, 59 FRONT_EFFICIENCY, 60 BACK_ABSORBED_FLUX, 61 BACK_INCOMING_FLUX, 62 BACK_ABSORBED_FIELD_GAIN, 63 BACK_ABSORBED_ATM_GAIN, 64 BACK_EFFICIENCY, 65 RECEIVER_RESULTS_COUNT__ 66 }; 67 68 enum primary_result_type { 69 FIRST_PRIMARY_RESULT, 70 PRIMARY_COS = FIRST_PRIMARY_RESULT, 71 PRIMARY_SHADOW, 72 PRIMARY_RESULTS_COUNT__ 73 }; 74 75 struct counts { 76 unsigned long global, receiver, primary, realisation, failed; 77 }; 78 79 static int 80 counts_ok(const struct counts* ref, const struct counts* c) 81 { 82 CHK(ref->global == GLOBAL_RESULTS_COUNT__); 83 CHK(c->global == GLOBAL_RESULTS_COUNT__); 84 return ref->receiver == c->receiver 85 && ref->primary == c->primary 86 && ref->failed >= c->failed; 87 } 88 89 #define MAX_LINE_LEN 2048 90 91 static const char 92 sundir_header [] = "#--- Sun direction:"; 93 94 #define IS_NEW_BLOCK(Line, Header) (!strncmp((Line), (Header), strlen(Header))) 95 96 97 #ifdef COMPILER_CL 98 /* mkstemp extracted from libc/sysdeps/posix/tempname.c. Copyright 99 * (C) 1991-1999, 2000, 2001, 2006 Free Software Foundation, Inc. 100 * 101 * The GNU C Library is free software; you can redistribute it and/or 102 * modify it under the terms of the GNU Lesser General Public 103 * License as published by the Free Software Foundation; either 104 * version 2.1 of the License, or (at your option) any later version. */ 105 106 static const char letters [] = 107 "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789"; 108 109 /* Generate a temporary file name based on TMPL. TMPL must match the 110 * rules for mk[s]temp (i.e. end in "XXXXXX"). The name constructed 111 * does not exist at the time of the call to mkstemp. TMPL is 112 * overwritten with the result. */ 113 int 114 mkstemp(char *tmpl) 115 { 116 size_t len; 117 char *XXXXXX; 118 static unsigned long long value; 119 unsigned long long random_time_bits; 120 unsigned int count; 121 int fd = -1; 122 int save_errno = errno; 123 124 /* A lower bound on the number of temporary files to attempt to 125 * generate. The maximum total number of temporary file names that 126 * can exist for a given template is 62**6. It should never be 127 * necessary to try all these combinations. Instead if a reasonable 128 * number of names is tried (we define reasonable as 62**3) fail to 129 * give the system administrator the chance to remove the problems. */ 130 #define ATTEMPTS_MIN (62 * 62 * 62) 131 132 /* The number of times to attempt to generate a temporary file. To 133 * conform to POSIX, this must be no smaller than TMP_MAX. */ 134 #if ATTEMPTS_MIN < TMP_MAX 135 unsigned int attempts = TMP_MAX; 136 #else 137 unsigned int attempts = ATTEMPTS_MIN; 138 #endif 139 140 len = strlen(tmpl); 141 if (len < 6 || strcmp(&tmpl[len - 6], "XXXXXX")) { 142 errno = EINVAL; 143 return -1; 144 } 145 146 /* This is where the Xs start. */ 147 XXXXXX = &tmpl[len - 6]; 148 149 /* Get some more or less random data. */ 150 { 151 SYSTEMTIME stNow; 152 FILETIME ftNow; 153 154 /* get system time */ 155 GetSystemTime(&stNow); 156 stNow.wMilliseconds = 500; 157 if (!SystemTimeToFileTime(&stNow, &ftNow)) { 158 errno = -1; 159 return -1; 160 } 161 162 random_time_bits = (((unsigned long long)ftNow.dwHighDateTime << 32) 163 | (unsigned long long)ftNow.dwLowDateTime); 164 } 165 value += random_time_bits ^ (unsigned long long)GetCurrentThreadId(); 166 167 for(count = 0; count < attempts; value += 7777, ++count) { 168 unsigned long long v = value; 169 170 /* Fill in the random bits. */ 171 XXXXXX[0] = letters[v % 62]; 172 v /= 62; 173 XXXXXX[1] = letters[v % 62]; 174 v /= 62; 175 XXXXXX[2] = letters[v % 62]; 176 v /= 62; 177 XXXXXX[3] = letters[v % 62]; 178 v /= 62; 179 XXXXXX[4] = letters[v % 62]; 180 v /= 62; 181 XXXXXX[5] = letters[v % 62]; 182 183 fd = open(tmpl, O_RDWR|O_CREAT|O_EXCL, S_IREAD|S_IWRITE); 184 if (fd >= 0) { 185 errno = save_errno; 186 return fd; 187 } 188 else if (errno != EEXIST) 189 return -1; 190 } 191 192 /* We got out of the loop because we ran out of combinations to try. */ 193 errno = EEXIST; 194 return -1; 195 } 196 #endif 197 198 static int 199 read_line(char* line, size_t max_line_len, FILE* stream) 200 { 201 ASSERT(stream && line && max_line_len); 202 line = fgets(line, (int)max_line_len, stream); 203 if(!line) return 0; 204 CHK(strlen(line) + 1 < max_line_len); 205 return 1; 206 } 207 208 static int 209 get_angles_and_counts 210 (FILE* file, 211 double angles[2], 212 struct counts* counts) 213 { 214 char line[MAX_LINE_LEN]; 215 int r; 216 217 CHK(file != NULL); 218 CHK(angles != NULL); 219 CHK(counts != NULL); 220 221 /* Get sun dir */ 222 r = read_line(line, sizeof(line), file); 223 if (!r) { 224 CHK(feof(file) == 1); 225 return 0; 226 } 227 CHK(IS_NEW_BLOCK(line, sundir_header) == 1); 228 CHK(sscanf(line+strlen(sundir_header), "%lg%lg", &angles[0], &angles[1]) == 2); 229 230 /* Get counts */ 231 CHK(read_line(line, sizeof(line), file) == 1); 232 CHK( 233 sscanf(line, 234 "%lu %lu %lu %lu %lu", 235 &counts->global, &counts->receiver, &counts->primary, 236 &counts->realisation, &counts->failed) == 5); 237 return 1; 238 } 239 240 static void 241 read_global(FILE* file, double* E, double* SE) 242 { 243 char line[MAX_LINE_LEN]; 244 CHK(read_line(line, sizeof(line), file) == 1); 245 CHK(sscanf(line, "%lg %lg", E, SE) == 2); 246 } 247 248 static void 249 read_recv(FILE* file, char name[], double E[], double SE[]) 250 { 251 char line[MAX_LINE_LEN]; 252 253 CHK(file != NULL); 254 CHK(name != NULL); 255 CHK(E != NULL); 256 CHK(SE != NULL); 257 258 CHK(read_line(line, sizeof(line), file) == 1); 259 CHK( 260 sscanf(line, 261 "%s %*u %*g " 262 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg " 263 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", 264 name, /* ID, area */ 265 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 266 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 267 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 268 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 269 &E[FRONT_EFFICIENCY], &SE[FRONT_EFFICIENCY], 270 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 271 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 272 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 273 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN], 274 &E[BACK_EFFICIENCY], &SE[BACK_EFFICIENCY]) == 275 2 * RECEIVER_RESULTS_COUNT__ + 1); 276 } 277 278 static void 279 read_primary 280 (FILE* file, char name[], double* area, double E[], double SE[]) 281 { 282 char line[MAX_LINE_LEN]; 283 284 CHK(file != NULL); 285 CHK(area != NULL); 286 CHK(E != NULL); 287 CHK(SE != NULL); 288 289 CHK(read_line(line, sizeof(line), file) == 1); 290 CHK( 291 sscanf(line, 292 "%s %*u " 293 "%lg %*u " 294 "%lg %lg %lg %lg\n", 295 name, /* ID */ 296 area, /* count, */ 297 &E[PRIMARY_COS], &SE[PRIMARY_COS], 298 &E[PRIMARY_SHADOW], &SE[PRIMARY_SHADOW]) == 299 2 * PRIMARY_RESULTS_COUNT__ + 2); 300 } 301 302 303 static void 304 read_recvXprim 305 (FILE* file, 306 unsigned long* rcv_id, 307 unsigned long* prim_id, 308 double E[], 309 double SE[]) 310 { 311 char line[MAX_LINE_LEN]; 312 313 CHK(file != NULL); 314 CHK(rcv_id != NULL); 315 CHK(prim_id != NULL); 316 CHK(E != NULL); 317 CHK(SE != NULL); 318 319 CHK(read_line(line, sizeof(line), file) == 1); 320 CHK( 321 sscanf(line, 322 "%lu %lu " 323 "%lg %lg %lg %lg %lg %lg %lg %lg " 324 "%lg %lg %lg %lg %lg %lg %lg %lg", 325 rcv_id, prim_id, 326 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 327 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 328 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 329 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 330 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 331 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 332 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 333 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN]) == 334 2 * (RECEIVER_RESULTS_COUNT__ - 2 /* efficiencies not read */) + 2); 335 } 336 337 static void 338 compute_estimate_intersection 339 (double intersection[2], 340 const double scale, 341 const double E0, 342 const double SE0, 343 const double E1, 344 const double SE1) 345 { 346 double interval0[2], interval1[2]; 347 CHK(scale > 0); 348 interval0[0] = E0 - scale*SE0; 349 interval0[1] = E0 + scale*SE0; 350 interval1[0] = E1 - scale*SE1; 351 interval1[1] = E1 + scale*SE1; 352 intersection[0] = MMAX(interval0[0], interval1[0]); 353 intersection[1] = MMIN(interval0[1], interval1[1]); 354 } 355 356 static void 357 check_estimate 358 (double ref_E, 359 double ref_SE, 360 double test_E, 361 double test_SE) 362 { 363 if(ref_E == -1) { 364 CHK(ref_SE == -1); 365 CHK(test_E == -1); 366 CHK(test_SE == -1); 367 } else { 368 double interval[2]; 369 CHK(ref_SE >= 0); 370 CHK(test_E >= 0); 371 CHK(test_SE >= 0); 372 if(!ref_SE) ref_SE = ref_E / 1000.0; 373 if(!test_SE) test_SE = test_E / 1000.0; 374 compute_estimate_intersection(interval, 2, ref_E, ref_SE, test_E, test_SE); 375 CHK(interval[0] <= interval[1]); 376 } 377 } 378 379 static void 380 check_1_reference 381 (FILE* ref_file, 382 FILE* test_file, 383 const struct counts* counts) 384 { 385 unsigned n; 386 387 CHK(ref_file != NULL); 388 CHK(test_file != NULL); 389 CHK(counts != NULL); 390 391 /* both files' pointer are just past the new bloc header */ 392 393 for(n = 0; n < counts->global; n++) { 394 double reference_E, reference_SE, test_E, test_SE; 395 read_global(ref_file, &reference_E, &reference_SE); 396 read_global(test_file, &test_E, &test_SE); 397 check_estimate(reference_E, reference_SE, test_E, test_SE); 398 } 399 for(n = 0; n < counts->receiver; n++) { 400 char ref_rcv_name[MAX_LINE_LEN], test_rcv_name[MAX_LINE_LEN]; 401 double reference_E[RECEIVER_RESULTS_COUNT__]; 402 double reference_SE[RECEIVER_RESULTS_COUNT__]; 403 double test_E[RECEIVER_RESULTS_COUNT__]; 404 double test_SE[RECEIVER_RESULTS_COUNT__]; 405 enum receiver_result_type r; 406 407 read_recv(ref_file, ref_rcv_name, reference_E, reference_SE); 408 read_recv(test_file, test_rcv_name, test_E, test_SE); 409 CHK(strcmp(ref_rcv_name, test_rcv_name) == 0); 410 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 411 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 412 } 413 } 414 for(n = 0; n < counts->primary; n++) { 415 char ref_prim_name[MAX_LINE_LEN], test_prim_name[MAX_LINE_LEN]; 416 double reference_E[PRIMARY_RESULTS_COUNT__]; 417 double reference_SE[PRIMARY_RESULTS_COUNT__]; 418 double test_E[PRIMARY_RESULTS_COUNT__]; 419 double test_SE[PRIMARY_RESULTS_COUNT__]; 420 double ref_area, test_area; 421 enum primary_result_type r; 422 423 read_primary(ref_file, ref_prim_name, &ref_area, reference_E, reference_SE); 424 read_primary(test_file, test_prim_name, &test_area, test_E, test_SE); 425 check_estimate(ref_area, 0, test_area, 0); 426 CHK(strcmp(ref_prim_name, test_prim_name) == 0); 427 FOR_EACH(r, FIRST_PRIMARY_RESULT, PRIMARY_RESULTS_COUNT__) { 428 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 429 } 430 } 431 for(n = 0; n < counts->receiver * counts->primary; n++) { 432 double reference_E[RECEIVER_RESULTS_COUNT__]; 433 double reference_SE[RECEIVER_RESULTS_COUNT__]; 434 double test_E[RECEIVER_RESULTS_COUNT__]; 435 double test_SE[RECEIVER_RESULTS_COUNT__]; 436 unsigned long ref_rcv_id, ref_prim_id; 437 unsigned long test_rcv_id, test_prim_id; 438 439 enum receiver_result_type r; 440 read_recvXprim(ref_file, &ref_rcv_id, &ref_prim_id, reference_E, reference_SE); 441 read_recvXprim(test_file, &test_rcv_id, &test_prim_id, test_E, test_SE); 442 /* we rely on the order of outputs */ 443 CHK(ref_rcv_id == test_rcv_id); 444 CHK(ref_prim_id == test_prim_id); 445 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 446 if (r == FRONT_EFFICIENCY || r == BACK_EFFICIENCY) 447 continue; /* not read */ 448 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 449 } 450 } 451 } 452 453 static FINLINE int 454 create_tmp_file(char* name, const size_t max_sizeof_name) 455 { 456 const char* template = "solstice_tmp_file_XXXXXX"; 457 int fd; 458 CHK(name != NULL); 459 CHK(strlen(template)+1 <= max_sizeof_name-1); 460 strcpy(name, template); 461 fd = mkstemp(name); 462 CHK(fd != -1); 463 return fd; 464 } 465 466 static void 467 do_check(const char* binary, const char* dir, const char* base_name) 468 { 469 char ref_file_name[128]; 470 FILE* ref_file; 471 struct counts ref_counts, test_counts; 472 int n; 473 int err; 474 ASSERT(base_name); 475 476 n = snprintf(ref_file_name, sizeof(ref_file_name), "%s%s.ref", dir, base_name); 477 CHK((size_t)n < sizeof(ref_file_name)); 478 479 ref_file = fopen(ref_file_name, "r"); 480 CHK(ref_file != NULL); 481 482 while(!feof(ref_file)) { 483 char cmd[512]; 484 char test_file_name[128]; 485 double ref_sun_angles[2], test_sun_angles[2]; 486 FILE* test_file = NULL; 487 int fd = -1; 488 489 if (!get_angles_and_counts(ref_file, ref_sun_angles, &ref_counts)) 490 break; /* EOF */ 491 492 fd = create_tmp_file(test_file_name, sizeof(test_file_name)); 493 test_file = fdopen(fd, "r"); 494 CHK(test_file != NULL); 495 496 n = snprintf(cmd, sizeof(cmd), 497 "%s -o %s -f -D %g,%g -n %lu -R %s%s_receiver.yaml %s%s.yaml", 498 binary, test_file_name, SPLIT2(ref_sun_angles), ref_counts.realisation, 499 dir, base_name, dir, base_name); 500 CHK((unsigned)n < sizeof(cmd)); 501 502 err = system(cmd); 503 CHK(err == 0); 504 505 get_angles_and_counts(test_file, test_sun_angles, &test_counts); 506 CHK(d2_eq(ref_sun_angles, test_sun_angles) == 1); 507 CHK(counts_ok(&ref_counts, &test_counts) == 1); 508 check_1_reference(ref_file, test_file, &ref_counts); 509 510 fclose(test_file); 511 remove(test_file_name); 512 } 513 } 514 515 int 516 main(int argc, char** argv) 517 { 518 int err = 0; 519 520 if(argc != 4) { 521 printf("Usage: %s <solstice-binary> <file-path> <file-base-name>\n", argv[0]); 522 goto error; 523 } 524 525 do_check(argv[1], argv[2], argv[3]); 526 527 exit: 528 return err; 529 error: 530 err = 1; 531 goto exit; 532 } 533