simu

Ice hockey final standings simulator
git clone https://git.inz.fi/simu
Log | Files | Refs

simu.c (12514B)


      1 /*
      2  * Input:
      3  *  <number of teams>
      4  *  Team1
      5  *  Team2
      6  *
      7  *  <number of played games>
      8  *  <index of home team> <index of away team> <home score> <away score> <OT/SO>
      9  *
     10  *  <number of remaining games>
     11  *  <index of home team> <index of away team>
     12  *
     13  * Output:
     14  *  Team1 <current games> <current points> <elo rating> <total games> <expected points> <probab finish 1st> <probab finish 2nd> ...
     15  *  Team2 <current games> <current points> <elo rating> <total games> <expected points> <probab finish 1st> <probab finish 2nd> ...
     16  */
     17 #define _POSIX_C_SOURCE 200811L
     18 #include <limits.h>
     19 #include <math.h>
     20 #include <pthread.h>
     21 #include <stdbool.h>
     22 #include <stdint.h>
     23 #include <stdio.h>
     24 #include <stdlib.h>
     25 #include <string.h>
     26 #include <unistd.h>
     27 
     28 #if defined(USE_LEHMER)
     29 #define R_MAX UINT32_MAX
     30 struct rand_t {
     31 	__uint128_t state;
     32 };
     33 
     34 static void getseed(struct rand_t *r)
     35 {
     36 	struct timespec s;
     37 	clock_gettime(CLOCK_MONOTONIC, &s);
     38 	r->state = s.tv_nsec;
     39 	r->state <<= 64;
     40 	r->state |= s.tv_sec | 1;
     41 }
     42 
     43 static inline uint32_t getrand(struct rand_t *r)
     44 {
     45 	r->state *= 0xda942042e4dd58b5;
     46 	return r->state >> 64;
     47 }
     48 
     49 #elif defined(USE_URANDOM)
     50 #define R_MAX UINT32_MAX
     51 struct rand_t {
     52 	int fd;
     53 	uint32_t rbuf[4096];
     54 	uint32_t *r;
     55 };
     56 
     57 static void getseed(struct rand_t *r)
     58 {
     59 	r->fd = open("/dev/urandom", O_RDONLY);
     60 	r->r = 1[&r->rbuf];
     61 }
     62 
     63 static inline uint32_t getrand(struct rand_t *r)
     64 {
     65 	if (r->r == 1[&r->rbuf]) {
     66 		read(r->fd, r->rbuf, sizeof(r->rbuf));
     67 		r->r = r->rbuf;
     68 	}
     69 	return *r->r++;
     70 }
     71 #else
     72 struct rand_t {
     73 	unsigned int seed;
     74 };
     75 #define R_MAX RAND_MAX
     76 static void getseed(struct rand_t *r)
     77 {
     78 	struct timespec s;
     79 	clock_gettime(CLOCK_MONOTONIC, &s);
     80 	r->seed = s.tv_sec ^ s.tv_nsec;
     81 }
     82 
     83 static inline uint32_t getrand(struct rand_t *r)
     84 {
     85 	return rand_r(&r->seed);
     86 }
     87 #endif
     88 
     89 static inline void urand(uint32_t limit, bool rv[static 2], struct rand_t *r)
     90 {
     91 	uint32_t v = getrand(r);
     92 	rv[0] = (uint32_t)v < limit;
     93 	rv[1] = (uint32_t)R_MAX - v < limit;
     94 }
     95 
     96 struct team {
     97 	size_t id;
     98 	char name[128];
     99 	unsigned points;
    100 	unsigned games;
    101 	unsigned wins;
    102 	double elo;
    103 	double elo0;
    104 	long long unsigned pointsum;
    105 	long long unsigned gamessum;
    106 	size_t *poscounts;
    107 };
    108 
    109 struct game {
    110 	size_t t[2];
    111 	uint32_t expected;
    112 };
    113 
    114 struct standing {
    115 	size_t ti;
    116 	unsigned points;
    117 	unsigned wins;
    118 	unsigned games;
    119 	unsigned random;
    120 };
    121 
    122 int pointscmp(const void *a, const void *b)
    123 {
    124 	int r;
    125 	const struct standing *at = a;
    126 	const struct standing *bt = b;
    127 
    128 	if ((r = bt->points * at->games - at->points * bt->games))
    129 		return r;
    130 	if ((r = bt->wins - at->wins))
    131 		return r;
    132 	return bt->random - at->random;
    133 }
    134 
    135 int teampointscmp(const void *a, const void *b)
    136 {
    137 	int r;
    138 	const struct team *at = a;
    139 	const struct team *bt = b;
    140 
    141 	if ((r = bt->points * at->games - at->points * bt->games))
    142 		return r;
    143 	if ((r = bt->wins - at->wins))
    144 		return r;
    145 	return 0;
    146 }
    147 
    148 int teamcmp(const void *a, const void *b)
    149 {
    150 	const struct team *at = a;
    151 	const struct team *bt = b;
    152 
    153 	if (at->pointsum * bt->gamessum > bt->pointsum * at->gamessum)
    154 		return -1;
    155 	if (bt->pointsum * at->gamessum > at->pointsum * bt->gamessum)
    156 		return 1;
    157 	return 0;
    158 }
    159 
    160 void isort(void *data, size_t n, size_t sz, int (*cmp)(const void *a, const void *b))
    161 {
    162 	char (*d)[sz][1] = data;
    163 	char tmp[sz];
    164 	size_t i, j;
    165 
    166 	for (i = 1; i < n; i++) {
    167 		if (cmp(d[i], d[i - 1]) >= 0)
    168 			continue;
    169 
    170 		memcpy(tmp, d[i], sizeof(tmp));
    171 
    172 		for (j = i - 1; j > 0 && cmp(tmp, d[j - 1]) < 0; j--);
    173 
    174 		memmove(d[j + 1], d[j], *d[i] - *d[j]);
    175 		memcpy(d[j], tmp, sizeof(tmp));
    176 	}
    177 }
    178 
    179 void simulate(struct team *teams,
    180 	      size_t n_teams,
    181 	      const struct game *games,
    182 	      size_t n_games,
    183 	      uint32_t otprob,
    184 	      size_t iterations,
    185 	      struct rand_t *seedp)
    186 {
    187 	struct standing s0[n_teams];
    188 	struct standing standings[2][n_teams];
    189 	size_t i;
    190 
    191 	for (i = 0; i < n_teams; i++) {
    192 		s0[i].ti = i;
    193 		s0[i].points = teams[i].points;
    194 		s0[i].wins = teams[i].wins;
    195 		s0[i].games = teams[i].games;
    196 	}
    197 
    198 	for (i = 0; i < iterations / 2; i++) {
    199 		size_t j;
    200 		size_t k;
    201 
    202 		memcpy(standings[0], s0, sizeof(*standings));
    203 		memcpy(standings[1], s0, sizeof(*standings));
    204 
    205 		for (j = 0; j < n_games; j++) {
    206 			bool hw[2];
    207 			bool ot[2];
    208 
    209 			urand(games[j].expected, hw, seedp);
    210 			urand(otprob, ot, seedp);
    211 
    212 			for (k = 0; k < 2; k++) {
    213 				size_t wi = games[j].t[!hw[k]];
    214 				size_t li = games[j].t[hw[k]];
    215 
    216 				standings[k][wi].points += 3 - ot[k];
    217 				standings[k][li].points += ot[k];
    218 				standings[k][wi].wins += !ot[k];
    219 				standings[k][wi].games++;
    220 				standings[k][li].games++;
    221 			}
    222 		}
    223 
    224 		for (k = 0; k < 2; k++) {
    225 			for (j = 0; j < n_teams; j++)
    226 				standings[k][j].random = getrand(seedp) % (INT_MAX / 2);
    227 			isort(standings[k], n_teams, sizeof(*standings[k]), pointscmp);
    228 
    229 			for (j = 0; j < n_teams; j++) {
    230 				teams[standings[k][j].ti].pointsum += standings[k][j].points;
    231 				teams[standings[k][j].ti].gamessum += standings[k][j].games;
    232 				teams[standings[k][j].ti].poscounts[j]++;
    233 			}
    234 		}
    235 	}
    236 
    237 }
    238 
    239 void simulate_winall(struct team *teams,
    240 		     size_t n_teams,
    241 		     const struct game *games,
    242 		     size_t n_games,
    243 		     uint32_t otprob,
    244 		     size_t iterations,
    245 		     size_t team_id,
    246 		     bool win,
    247 		     struct rand_t *seedp)
    248 {
    249 	struct team teams_cpy[n_teams];
    250 	struct game games_cpy[n_games];
    251 
    252 	size_t r;
    253 	size_t w;
    254 
    255 	memcpy(teams_cpy, teams, sizeof(teams_cpy));
    256 
    257 	for (w = r = 0; r < n_games; r++) {
    258 		if (games[r].t[0] == team_id || games[r].t[1] == team_id) {
    259 			size_t opp = games[r].t[0] + games[r].t[1] - team_id;
    260 			teams_cpy[team_id].points += win * 3;
    261 			teams_cpy[team_id].games++;
    262 			teams_cpy[team_id].wins += win;
    263 
    264 			teams_cpy[opp].points += !win * 3;
    265 			teams_cpy[opp].games++;
    266 			teams_cpy[opp].wins += !win;
    267 		} else {
    268 			games_cpy[w++] = games[r];
    269 		}
    270 	}
    271 
    272 	simulate(teams_cpy, n_teams, games_cpy, w, otprob, iterations, seedp);
    273 }
    274 
    275 void simulate_winloseall(struct team *teams,
    276 			 size_t n_teams,
    277 			 size_t team,
    278 			 const struct game *games,
    279 			 size_t n_games,
    280 			 uint32_t otprob,
    281 			 size_t iterations,
    282 			 struct rand_t *seedp)
    283 {
    284 	struct team teams_cpy[n_teams];
    285 	size_t poscounts[n_teams][n_teams];
    286 	size_t e;
    287 	size_t i;
    288 
    289 	memcpy(teams_cpy, teams, sizeof(teams_cpy));
    290 	memset(poscounts, 0, sizeof(poscounts));
    291 
    292 	for (i = 0; i < n_teams; i++)
    293 		teams_cpy[i].poscounts = poscounts[i];
    294 
    295 	simulate_winall(teams_cpy, n_teams,
    296 			games, n_games,
    297 			otprob,
    298 			iterations,
    299 			team, true,
    300 			seedp);
    301 	simulate_winall(teams_cpy, n_teams,
    302 			games, n_games,
    303 			otprob,
    304 			iterations,
    305 			team, false,
    306 			seedp);
    307 
    308 	for (i = 0; i < n_teams; i++)
    309 		if (teams_cpy[team].poscounts[i])
    310 			break;
    311 	for (e = n_teams; e > 0; e--)
    312 		if (teams_cpy[team].poscounts[e - 1])
    313 			break;
    314 	for (; i < e; i++)
    315 		teams[team].poscounts[i] = 1;
    316 }
    317 
    318 struct thread_data {
    319 	struct team *teams;
    320 	bool *winloseall;
    321 	size_t n_teams;
    322 
    323 	struct game *games;
    324 	size_t n_games;
    325 
    326 	uint32_t otprob;
    327 
    328 	size_t iterations;
    329 	size_t iterated;
    330 	pthread_mutex_t mutex;
    331 };
    332 
    333 void *simulate_thread(void *data)
    334 {
    335 	struct thread_data *td = data;
    336 	struct team teams_cpy[td->n_teams];
    337 	size_t i;
    338 	size_t iters;
    339 	size_t poscounts[td->n_teams][td->n_teams];
    340 	size_t n = 0;
    341 	struct rand_t seedp;
    342 
    343 	getseed(&seedp);
    344 
    345 	pthread_mutex_lock(&td->mutex);
    346 	memcpy(teams_cpy, td->teams, sizeof(teams_cpy));
    347 	pthread_mutex_unlock(&td->mutex);
    348 
    349 	memset(poscounts, 0, sizeof(poscounts));
    350 
    351 	for (i = 0; i < td->n_teams; i++) {
    352 		teams_cpy[i].poscounts = poscounts[i];
    353 		pthread_mutex_lock(&td->mutex);
    354 		if (!td->winloseall[i]) {
    355 			td->winloseall[i] = true;
    356 			pthread_mutex_unlock(&td->mutex);
    357 
    358 			simulate_winloseall(teams_cpy, td->n_teams, i,
    359 					    td->games, td->n_games,
    360 					    td->otprob,
    361 					    td->iterations / td->n_teams / 2,
    362 					    &seedp);
    363 		} else {
    364 			pthread_mutex_unlock(&td->mutex);
    365 		}
    366 	}
    367 
    368 	for (;;) {
    369 		pthread_mutex_lock(&td->mutex);
    370 
    371 		if (td->iterations - td->iterated > 100000)
    372 			iters = 100000;
    373 		else
    374 			iters = td->iterations - td->iterated;
    375 
    376 		td->iterated += iters;
    377 		pthread_mutex_unlock(&td->mutex);
    378 
    379 		if (!iters)
    380 			break;
    381 		n += iters;
    382 
    383 		simulate(teams_cpy, td->n_teams,
    384 			 td->games, td->n_games,
    385 			 td->otprob,
    386 			 iters, &seedp);
    387 	}
    388 
    389 	pthread_mutex_lock(&td->mutex);
    390 	for (i = 0; i < td->n_teams; i++) {
    391 		size_t j;
    392 		td->teams[i].pointsum += teams_cpy[i].pointsum;
    393 		td->teams[i].gamessum += teams_cpy[i].gamessum;
    394 
    395 		for (j = 0; j < td->n_teams; j++)
    396 			td->teams[i].poscounts[j] += teams_cpy[i].poscounts[j];
    397 	}
    398 	pthread_mutex_unlock(&td->mutex);
    399 
    400 	return (void *)n;
    401 }
    402 
    403 const char *ltrim(const char *s)
    404 {
    405 	while (*s++ == ' ');
    406 	return s - 1;
    407 }
    408 
    409 int main(int argc, char **argv)
    410 {
    411 	int opt;
    412 	double homeadv = 58;
    413 	double elok = 32;
    414 	uint32_t otprob = 0.23 * R_MAX;
    415 	size_t n;
    416 	size_t m;
    417 	size_t i;
    418 	size_t j;
    419 	size_t iterations = 1000000;
    420 	size_t threads = 1;
    421 
    422 	while ((opt = getopt(argc, argv, "i:t:a:o:")) != -1) {
    423 		switch (opt) {
    424 		case 'i':
    425 			iterations = strtoull(optarg, NULL, 10);
    426 			break;
    427 
    428 		case 't':
    429 			threads = strtoull(optarg, NULL, 10);
    430 			break;
    431 
    432 		case 'a':
    433 			homeadv = strtod(optarg, NULL);
    434 			break;
    435 
    436 		case 'o':
    437 			otprob = strtod(optarg, NULL) * R_MAX;
    438 			break;
    439 
    440 		case 'k':
    441 			elok = strtod(optarg, NULL);
    442 			break;
    443 
    444 		default:
    445 			exit(1);
    446 			break;
    447 		}
    448 	}
    449 
    450 	iterations += iterations & 1;
    451 
    452 	if (scanf(" %zu", &n) < 1)
    453 		return 1;
    454 
    455 	struct team teams[n];
    456 	size_t teammap[n];
    457 	size_t poscounts[n][n];
    458 	memset(poscounts, 0, sizeof(poscounts));
    459 
    460 	for (i = 0; i < n; i++) {
    461 		if (scanf(" %127s", teams[i].name) < 1)
    462 			return 1;
    463 		teams[i].id = i;
    464 		teams[i].points = 0;
    465 		teams[i].games = 0;
    466 		teams[i].wins = 0;
    467 		teams[i].pointsum = 0;
    468 		teams[i].gamessum = 0;
    469 		teams[i].poscounts = poscounts[i];
    470 		teams[i].elo = 2000;
    471 	}
    472 
    473 	if (scanf(" %zu", &m) < 1)
    474 		return 1;
    475 
    476 	for (i = 0; i < m; i++) {
    477 		size_t hi, ai;
    478 		unsigned hs, as;
    479 		bool ot = false;
    480 		char det[16] = "";
    481 
    482 		switch (scanf(" %zu %zu %u %u%15[^\n]", &hi, &ai, &hs, &as, det)) {
    483 		case 5:
    484 			ot = !strcmp(ltrim(det), "OT") ||
    485 			     !strcmp(ltrim(det), "SO");
    486 			break;
    487 
    488 		case 4:
    489 			break;
    490 
    491 		default:
    492 			exit(1);
    493 		}
    494 
    495 		double expected = 1. / (1 + pow(10, (teams[ai].elo - teams[hi].elo - homeadv) / 400.0));
    496 		double actual;
    497 
    498 		teams[hi].games++;
    499 		teams[ai].games++;
    500 
    501 		if (hs > as) {
    502 			if (!ot) {
    503 				actual = 1;
    504 				teams[hi].wins++;
    505 				teams[hi].points += 3;
    506 			} else {
    507 				actual = 2./3;
    508 				teams[hi].points += 2;
    509 				teams[ai].points += 1;
    510 			}
    511 		} else {
    512 			if (ot) {
    513 				actual = 1./3;
    514 				teams[hi].points += 1;
    515 				teams[ai].points += 2;
    516 			} else {
    517 				teams[ai].points += 3;
    518 				teams[ai].wins++;
    519 				actual = 0;
    520 			}
    521 		}
    522 
    523 		teams[hi].elo += elok * (actual - expected);
    524 		teams[ai].elo -= elok * (actual - expected);
    525 	}
    526 
    527 	isort(teams, 1[&teams] - teams, sizeof(*teams), teampointscmp);
    528 
    529 	for (i = 0; i < n; i++) {
    530 		teammap[teams[i].id] = i;
    531 		teams[i].elo0 = teams[i].elo;
    532 	}
    533 
    534 	if (scanf(" %zu", &m) < 1)
    535 		return 1;
    536 
    537 	struct game games[m];
    538 
    539 	for (i = 0; i < m; i++) {
    540 		size_t h, a;
    541 		if (scanf(" %zu %zu", &h, &a) < 2)
    542 			return 1;
    543 		games[i].t[0] = teammap[h];
    544 		games[i].t[1] = teammap[a];
    545 		games[i].expected = R_MAX / (1 + pow(10.0, (teams[games[i].t[1]].elo - (teams[games[i].t[0]].elo + homeadv)) / 400.0));
    546 
    547 		teams[games[i].t[0]].elo += (2000 - teams[games[i].t[0]].elo) * elok / 200;
    548 		teams[games[i].t[1]].elo += (2000 - teams[games[i].t[1]].elo) * elok / 200;
    549 	}
    550 
    551 	if (threads > 1) {
    552 		pthread_t pthrds[threads];
    553 		bool winloseall[n];
    554 		memset(winloseall, 0, sizeof(winloseall));
    555 		struct thread_data td = {
    556 			teams, winloseall, n,
    557 			games, m,
    558 			otprob,
    559 			iterations, 0,
    560 			PTHREAD_MUTEX_INITIALIZER
    561 		};
    562 		for (i = 0; i < threads; i++)
    563 			pthread_create(&pthrds[i], NULL, simulate_thread, &td);
    564 		iterations = 0;
    565 		for (i = 0; i < threads; i++) {
    566 			void *a;
    567 			pthread_join(pthrds[i], &a);
    568 			iterations += (size_t)a;
    569 		}
    570 	} else {
    571 		struct rand_t seedp;
    572 		getseed(&seedp);
    573 
    574 		for (i = 0; i < n; i++) {
    575 			simulate_winloseall(teams, n, i,
    576 					    games, m,
    577 					    otprob,
    578 					    iterations / n / 2,
    579 					    &seedp);
    580 		}
    581 
    582 		simulate(teams, n,
    583 			 games, m,
    584 			 otprob,
    585 			 iterations,
    586 			 &seedp);
    587 	}
    588 
    589 	isort(teams, n, sizeof(*teams), teamcmp);
    590 	for (i = 0; i < n; i++) {
    591 		printf("%s %u %u %llu %.1f %d", teams[i].name,
    592 		       teams[i].games, teams[i].points,
    593 		       teams[i].gamessum / iterations,
    594 		       (double)teams[i].pointsum / iterations,
    595 		       (int)round(teams[i].elo0));
    596 		for (j = 0; j < n; j++) {
    597 			if (teams[i].poscounts[j] == iterations + 1)
    598 				printf(" +");
    599 			else if (teams[i].poscounts[j])
    600 				printf(" %e", (teams[i].poscounts[j] - 1) * 100.0 / iterations);
    601 			else
    602 				printf(" -");
    603 		}
    604 		printf("\n");
    605 	}
    606 
    607 	return 1;
    608 }