test_solstice_simulation.c (11672B)
1 /* Copyright (C) 2018-2026 |Méso|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 <stdio.h> 24 #include <stdlib.h> 25 #include <string.h> 26 27 enum side { 28 FRONT, 29 BACK 30 }; 31 32 enum global_result_type { 33 GLOBAL_POTENTIAL, 34 GLOBAL_ABSORBED, 35 GLOBAL_COS, 36 GLOBAL_SHADOW, 37 GLOBAL_MISSING, 38 GLOBAL_ATMOSPHERE, 39 GLOBAL_REFLECTIVITY, 40 GLOBAL_RESULTS_COUNT__ 41 }; 42 43 enum receiver_result_type { 44 FIRST_RECEIVER_RESULT, 45 FRONT_ABSORBED_FLUX = FIRST_RECEIVER_RESULT, 46 FRONT_INCOMING_FLUX, 47 FRONT_ABSORBED_FIELD_GAIN, 48 FRONT_ABSORBED_ATM_GAIN, 49 FRONT_EFFICIENCY, 50 BACK_ABSORBED_FLUX, 51 BACK_INCOMING_FLUX, 52 BACK_ABSORBED_FIELD_GAIN, 53 BACK_ABSORBED_ATM_GAIN, 54 BACK_EFFICIENCY, 55 RECEIVER_RESULTS_COUNT__ 56 }; 57 58 enum primary_result_type { 59 FIRST_PRIMARY_RESULT, 60 PRIMARY_COS = FIRST_PRIMARY_RESULT, 61 PRIMARY_SHADOW, 62 PRIMARY_RESULTS_COUNT__ 63 }; 64 65 struct counts { 66 unsigned long global, receiver, primary, realisation, failed; 67 }; 68 69 static int 70 counts_ok(const struct counts* ref, const struct counts* c) 71 { 72 CHK(ref->global == GLOBAL_RESULTS_COUNT__); 73 CHK(c->global == GLOBAL_RESULTS_COUNT__); 74 return ref->receiver == c->receiver 75 && ref->primary == c->primary 76 && ref->failed >= c->failed; 77 } 78 79 #define MAX_LINE_LEN 2048 80 81 static const char 82 sundir_header [] = "#--- Sun direction:"; 83 84 #define IS_NEW_BLOCK(Line, Header) (!strncmp((Line), (Header), strlen(Header))) 85 86 static int 87 read_line(char* line, size_t max_line_len, FILE* stream) 88 { 89 ASSERT(stream && line && max_line_len); 90 line = fgets(line, (int)max_line_len, stream); 91 if(!line) return 0; 92 CHK(strlen(line) + 1 < max_line_len); 93 return 1; 94 } 95 96 static int 97 get_angles_and_counts 98 (FILE* file, 99 double angles[2], 100 struct counts* counts) 101 { 102 char line[MAX_LINE_LEN]; 103 int r; 104 105 CHK(file != NULL); 106 CHK(angles != NULL); 107 CHK(counts != NULL); 108 109 /* Get sun dir */ 110 r = read_line(line, sizeof(line), file); 111 if (!r) { 112 CHK(feof(file) == 1); 113 return 0; 114 } 115 CHK(IS_NEW_BLOCK(line, sundir_header) == 1); 116 CHK(sscanf(line+strlen(sundir_header), "%lg%lg", &angles[0], &angles[1]) == 2); 117 118 /* Get counts */ 119 CHK(read_line(line, sizeof(line), file) == 1); 120 CHK( 121 sscanf(line, 122 "%lu %lu %lu %lu %lu", 123 &counts->global, &counts->receiver, &counts->primary, 124 &counts->realisation, &counts->failed) == 5); 125 return 1; 126 } 127 128 static void 129 read_global(FILE* file, double* E, double* SE) 130 { 131 char line[MAX_LINE_LEN]; 132 CHK(read_line(line, sizeof(line), file) == 1); 133 CHK(sscanf(line, "%lg %lg", E, SE) == 2); 134 } 135 136 static void 137 read_recv(FILE* file, char name[], double E[], double SE[]) 138 { 139 char line[MAX_LINE_LEN]; 140 141 CHK(file != NULL); 142 CHK(name != NULL); 143 CHK(E != NULL); 144 CHK(SE != NULL); 145 146 CHK(read_line(line, sizeof(line), file) == 1); 147 CHK( 148 sscanf(line, 149 "%s %*u %*g " 150 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg " 151 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", 152 name, /* ID, area */ 153 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 154 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 155 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 156 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 157 &E[FRONT_EFFICIENCY], &SE[FRONT_EFFICIENCY], 158 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 159 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 160 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 161 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN], 162 &E[BACK_EFFICIENCY], &SE[BACK_EFFICIENCY]) == 163 2 * RECEIVER_RESULTS_COUNT__ + 1); 164 } 165 166 static void 167 read_primary 168 (FILE* file, char name[], double* area, double E[], double SE[]) 169 { 170 char line[MAX_LINE_LEN]; 171 172 CHK(file != NULL); 173 CHK(area != NULL); 174 CHK(E != NULL); 175 CHK(SE != NULL); 176 177 CHK(read_line(line, sizeof(line), file) == 1); 178 CHK( 179 sscanf(line, 180 "%s %*u " 181 "%lg %*u " 182 "%lg %lg %lg %lg\n", 183 name, /* ID */ 184 area, /* count, */ 185 &E[PRIMARY_COS], &SE[PRIMARY_COS], 186 &E[PRIMARY_SHADOW], &SE[PRIMARY_SHADOW]) == 187 2 * PRIMARY_RESULTS_COUNT__ + 2); 188 } 189 190 191 static void 192 read_recvXprim 193 (FILE* file, 194 unsigned long* rcv_id, 195 unsigned long* prim_id, 196 double E[], 197 double SE[]) 198 { 199 char line[MAX_LINE_LEN]; 200 201 CHK(file != NULL); 202 CHK(rcv_id != NULL); 203 CHK(prim_id != NULL); 204 CHK(E != NULL); 205 CHK(SE != NULL); 206 207 CHK(read_line(line, sizeof(line), file) == 1); 208 CHK( 209 sscanf(line, 210 "%lu %lu " 211 "%lg %lg %lg %lg %lg %lg %lg %lg " 212 "%lg %lg %lg %lg %lg %lg %lg %lg", 213 rcv_id, prim_id, 214 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 215 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 216 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 217 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 218 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 219 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 220 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 221 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN]) == 222 2 * (RECEIVER_RESULTS_COUNT__ - 2 /* efficiencies not read */) + 2); 223 } 224 225 static void 226 compute_estimate_intersection 227 (double intersection[2], 228 const double scale, 229 const double E0, 230 const double SE0, 231 const double E1, 232 const double SE1) 233 { 234 double interval0[2], interval1[2]; 235 CHK(scale > 0); 236 interval0[0] = E0 - scale*SE0; 237 interval0[1] = E0 + scale*SE0; 238 interval1[0] = E1 - scale*SE1; 239 interval1[1] = E1 + scale*SE1; 240 intersection[0] = MMAX(interval0[0], interval1[0]); 241 intersection[1] = MMIN(interval0[1], interval1[1]); 242 } 243 244 static void 245 check_estimate 246 (double ref_E, 247 double ref_SE, 248 double test_E, 249 double test_SE) 250 { 251 if(ref_E == -1) { 252 CHK(ref_SE == -1); 253 CHK(test_E == -1); 254 CHK(test_SE == -1); 255 } else { 256 double interval[2]; 257 CHK(ref_SE >= 0); 258 CHK(test_E >= 0); 259 CHK(test_SE >= 0); 260 if(!ref_SE) ref_SE = ref_E / 1000.0; 261 if(!test_SE) test_SE = test_E / 1000.0; 262 compute_estimate_intersection(interval, 2, ref_E, ref_SE, test_E, test_SE); 263 CHK(interval[0] <= interval[1]); 264 } 265 } 266 267 static void 268 check_1_reference 269 (FILE* ref_file, 270 FILE* test_file, 271 const struct counts* counts) 272 { 273 unsigned n; 274 275 CHK(ref_file != NULL); 276 CHK(test_file != NULL); 277 CHK(counts != NULL); 278 279 /* both files' pointer are just past the new bloc header */ 280 281 for(n = 0; n < counts->global; n++) { 282 double reference_E, reference_SE, test_E, test_SE; 283 read_global(ref_file, &reference_E, &reference_SE); 284 read_global(test_file, &test_E, &test_SE); 285 check_estimate(reference_E, reference_SE, test_E, test_SE); 286 } 287 for(n = 0; n < counts->receiver; n++) { 288 char ref_rcv_name[MAX_LINE_LEN], test_rcv_name[MAX_LINE_LEN]; 289 double reference_E[RECEIVER_RESULTS_COUNT__]; 290 double reference_SE[RECEIVER_RESULTS_COUNT__]; 291 double test_E[RECEIVER_RESULTS_COUNT__]; 292 double test_SE[RECEIVER_RESULTS_COUNT__]; 293 enum receiver_result_type r; 294 295 read_recv(ref_file, ref_rcv_name, reference_E, reference_SE); 296 read_recv(test_file, test_rcv_name, test_E, test_SE); 297 CHK(strcmp(ref_rcv_name, test_rcv_name) == 0); 298 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 299 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 300 } 301 } 302 for(n = 0; n < counts->primary; n++) { 303 char ref_prim_name[MAX_LINE_LEN], test_prim_name[MAX_LINE_LEN]; 304 double reference_E[PRIMARY_RESULTS_COUNT__]; 305 double reference_SE[PRIMARY_RESULTS_COUNT__]; 306 double test_E[PRIMARY_RESULTS_COUNT__]; 307 double test_SE[PRIMARY_RESULTS_COUNT__]; 308 double ref_area, test_area; 309 enum primary_result_type r; 310 311 read_primary(ref_file, ref_prim_name, &ref_area, reference_E, reference_SE); 312 read_primary(test_file, test_prim_name, &test_area, test_E, test_SE); 313 check_estimate(ref_area, 0, test_area, 0); 314 CHK(strcmp(ref_prim_name, test_prim_name) == 0); 315 FOR_EACH(r, FIRST_PRIMARY_RESULT, PRIMARY_RESULTS_COUNT__) { 316 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 317 } 318 } 319 for(n = 0; n < counts->receiver * counts->primary; n++) { 320 double reference_E[RECEIVER_RESULTS_COUNT__]; 321 double reference_SE[RECEIVER_RESULTS_COUNT__]; 322 double test_E[RECEIVER_RESULTS_COUNT__]; 323 double test_SE[RECEIVER_RESULTS_COUNT__]; 324 unsigned long ref_rcv_id, ref_prim_id; 325 unsigned long test_rcv_id, test_prim_id; 326 327 enum receiver_result_type r; 328 read_recvXprim(ref_file, &ref_rcv_id, &ref_prim_id, reference_E, reference_SE); 329 read_recvXprim(test_file, &test_rcv_id, &test_prim_id, test_E, test_SE); 330 /* we rely on the order of outputs */ 331 CHK(ref_rcv_id == test_rcv_id); 332 CHK(ref_prim_id == test_prim_id); 333 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 334 if (r == FRONT_EFFICIENCY || r == BACK_EFFICIENCY) 335 continue; /* not read */ 336 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 337 } 338 } 339 } 340 341 static FINLINE int 342 create_tmp_file(char* name, const size_t max_sizeof_name) 343 { 344 const char* template = "solstice_tmp_file_XXXXXX"; 345 int fd; 346 CHK(name != NULL); 347 CHK(strlen(template)+1 <= max_sizeof_name-1); 348 strcpy(name, template); 349 fd = mkstemp(name); 350 CHK(fd != -1); 351 return fd; 352 } 353 354 static void 355 do_check(const char* binary, const char* dir, const char* base_name) 356 { 357 char ref_file_name[128]; 358 FILE* ref_file; 359 struct counts ref_counts, test_counts; 360 int n; 361 int err; 362 ASSERT(base_name); 363 364 n = snprintf(ref_file_name, sizeof(ref_file_name), "%s%s.ref", dir, base_name); 365 CHK((size_t)n < sizeof(ref_file_name)); 366 367 ref_file = fopen(ref_file_name, "r"); 368 CHK(ref_file != NULL); 369 370 while(!feof(ref_file)) { 371 char cmd[512]; 372 char test_file_name[128]; 373 double ref_sun_angles[2], test_sun_angles[2]; 374 FILE* test_file = NULL; 375 int fd = -1; 376 377 if (!get_angles_and_counts(ref_file, ref_sun_angles, &ref_counts)) 378 break; /* EOF */ 379 380 fd = create_tmp_file(test_file_name, sizeof(test_file_name)); 381 test_file = fdopen(fd, "r"); 382 CHK(test_file != NULL); 383 384 n = snprintf(cmd, sizeof(cmd), 385 "%s -o %s -f -D %g,%g -n %lu -R %s%s_receiver.yaml %s%s.yaml", 386 binary, test_file_name, SPLIT2(ref_sun_angles), ref_counts.realisation, 387 dir, base_name, dir, base_name); 388 CHK((unsigned)n < sizeof(cmd)); 389 390 err = system(cmd); 391 CHK(err == 0); 392 393 get_angles_and_counts(test_file, test_sun_angles, &test_counts); 394 CHK(d2_eq(ref_sun_angles, test_sun_angles) == 1); 395 CHK(counts_ok(&ref_counts, &test_counts) == 1); 396 check_1_reference(ref_file, test_file, &ref_counts); 397 398 fclose(test_file); 399 remove(test_file_name); 400 } 401 } 402 403 int 404 main(int argc, char** argv) 405 { 406 int err = 0; 407 408 if(argc != 4) { 409 printf("Usage: %s <solstice-binary> <file-path> <file-base-name>\n", argv[0]); 410 goto error; 411 } 412 413 do_check(argv[1], argv[2], argv[3]); 414 415 exit: 416 return err; 417 error: 418 err = 1; 419 goto exit; 420 } 421