simu

commit 8b28e56a62290385cd5605112d6d5bc9c6336167
Author: Santtu Lakkala <inz@inz.fi>
Date:   Thu,  7 Mar 2024 17:52:39 +0200

Initial import

Diffstat:
A
Makefile
 | 
8
++++++++
A
parsedata
 | 
31
+++++++++++++++++++++++++++++++
A
simu.c
 | 
621
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A
simuformat
 | 
92
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

4 files changed, 752 insertions(+), 0 deletions(-)
diff --git a/Makefile b/Makefile
@@ -0,0 +1,8 @@
+LDFLAGS = -lm
+CFLAGS = -W -Wall -std=c99 -pthread -DUSE_LEHMER
+TARGET = simu
+
+all: $(TARGET)
+
+$(TARGET): simu.c
+	$(CC) $(CFLAGS) -o $@ $< $(LDFLAGS)
diff --git a/parsedata b/parsedata
@@ -0,0 +1,31 @@
+#!/usr/bin/env perl
+
+use JSON;
+
+my $data = JSON->new->utf8->decode(join '', <<>>) || die "$!";
+my %teams;
+my @played;
+my @upcoming;
+
+for (@{$data}) {
+	my $ht = $_->{'homeTeam'}{'teamName'};
+	my $at = $_->{'awayTeam'}{'teamName'};
+	my $hid = $teams{$ht} // ($teams{$ht} = scalar %teams);
+	my $aid = $teams{$at} // ($teams{$at} = scalar %teams);
+
+	if ($_->{ended}) {
+		push @played, [$hid, $aid, $_->{'homeTeam'}{'goals'}, $_->{'awayTeam'}{'goals'}, $_->{'gameTime'}];
+	} else {
+		push @upcoming, [$hid, $aid];
+	}
+}
+
+printf("%d\n", scalar %teams);
+print join("\n", sort { $teams{$a} <=> $teams{$b} } keys %teams), "\n\n";
+
+printf("%d\n", scalar @played);
+print map { sprintf("%d %d %d %d %s\n", @$_[0..3], $_->[4] > 3600 ? $_->[4] == 3900 ? "SO" : "OT" : "") } @played;
+print "\n";
+
+printf("%d\n", scalar @upcoming);
+print map { (join ' ', @$_) . "\n" } @upcoming;
diff --git a/simu.c b/simu.c
@@ -0,0 +1,621 @@
+/*
+ * Input:
+ *  <number of teams>
+ *  Team1
+ *  Team2
+ *
+ *  <number of played games>
+ *  <index of home team> <index of away team> <home score> <away score> <OT/SO>
+ *
+ *  <number of remaining games>
+ *  <index of home team> <index of away team>
+ *
+ * Output:
+ *  Team1 <current games> <current points> <elo rating> <total games> <expected points> <probab finish 1st> <probab finish 2nd> ...
+ *  Team2 <current games> <current points> <elo rating> <total games> <expected points> <probab finish 1st> <probab finish 2nd> ...
+ */
+#define _POSIX_C_SOURCE 200811L
+#include <limits.h>
+#include <math.h>
+#include <pthread.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+#if defined(USE_LEHMER)
+#define R_MAX UINT32_MAX
+struct rand_t {
+	__uint128_t state;
+};
+
+static void getseed(struct rand_t *r)
+{
+	struct timespec s;
+	clock_gettime(CLOCK_MONOTONIC, &s);
+	r->state = s.tv_nsec;
+	r->state <<= 64;
+	r->state |= s.tv_sec | 1;
+}
+
+static inline uint32_t getrand(struct rand_t *r)
+{
+	r->state *= 0xda942042e4dd58b5;
+	return r->state >> 64;
+}
+
+#elif defined(USE_URANDOM)
+#define R_MAX UINT32_MAX
+struct rand_t {
+	int fd;
+	uint32_t rbuf[4096];
+	uint32_t *r;
+};
+
+static void getseed(struct rand_t *r)
+{
+	r->fd = open("/dev/urandom", O_RDONLY);
+	r->r = 1[&r->rbuf];
+}
+
+static inline uint32_t getrand(struct rand_t *r)
+{
+	if (r->r == 1[&r->rbuf]) {
+		read(r->fd, r->rbuf, sizeof(r->rbuf));
+		r->r = r->rbuf;
+	}
+	return *r->r++;
+}
+#else
+struct rand_t {
+	unsigned int seed;
+};
+#define R_MAX RAND_MAX
+static void getseed(struct rand_t *r)
+{
+	struct timespec s;
+	clock_gettime(CLOCK_MONOTONIC, &s);
+	r->seed = s.tv_sec ^ s.tv_nsec;
+}
+
+static inline uint32_t getrand(struct rand_t *r)
+{
+	return rand_r(&r->seed);
+}
+#endif
+
+static inline bool urand(uint32_t limit, bool *anti, struct rand_t *r)
+{
+	uint32_t v = getrand(r);
+	if (anti)
+		*anti = (uint32_t)R_MAX - v < limit;
+	return (uint32_t)v < limit;
+}
+
+struct team {
+	size_t id;
+	char name[128];
+	unsigned points;
+	unsigned games;
+	unsigned wins;
+	double elo;
+	long long unsigned pointsum;
+	long long unsigned gamessum;
+	size_t *poscounts;
+};
+
+struct game {
+	size_t t[2];
+	uint32_t expected;
+};
+
+struct standing {
+	size_t ti;
+	unsigned points;
+	unsigned wins;
+	unsigned games;
+	unsigned random;
+};
+
+int pointscmp(const void *a, const void *b)
+{
+	int r;
+	const struct standing *at = a;
+	const struct standing *bt = b;
+
+	if ((r = bt->points * at->games - at->points * bt->games))
+		return r;
+	if ((r = bt->wins - at->wins))
+		return r;
+	return bt->random - at->random;
+}
+
+int teampointscmp(const void *a, const void *b)
+{
+	int r;
+	const struct team *at = a;
+	const struct team *bt = b;
+
+	if ((r = bt->points * at->games - at->points * bt->games))
+		return r;
+	if ((r = bt->wins - at->wins))
+		return r;
+	return 0;
+}
+
+int teamcmp(const void *a, const void *b)
+{
+	const struct team *at = a;
+	const struct team *bt = b;
+
+	if (at->pointsum * bt->gamessum > bt->pointsum * at->gamessum)
+		return -1;
+	if (bt->pointsum * at->gamessum > at->pointsum * bt->gamessum)
+		return 1;
+	return 0;
+}
+
+void isort(void *data, size_t n, size_t sz, int (*cmp)(const void *a, const void *b))
+{
+	char (*d)[sz][1] = data;
+	char tmp[sz];
+	size_t i, j;
+
+	for (i = 1; i < n; i++) {
+		if (cmp(d[i], d[i - 1]) >= 0)
+			continue;
+
+		memcpy(tmp, d[i], sizeof(tmp));
+
+		for (j = i - 1; j > 0 && cmp(tmp, d[j - 1]) < 0; j--);
+
+		memmove(d[j + 1], d[j], *d[i] - *d[j]);
+		memcpy(d[j], tmp, sizeof(tmp));
+	}
+}
+
+void simulate(struct team *teams,
+	      size_t n_teams,
+	      const struct game *games,
+	      size_t n_games,
+	      uint32_t otprob,
+	      size_t iterations,
+	      struct rand_t *seedp)
+{
+	struct standing s0[n_teams];
+	struct standing standings[n_teams];
+	struct standing antistandings[n_teams];
+	size_t i;
+
+	for (i = 0; i < n_teams; i++) {
+		s0[i].ti = i;
+		s0[i].points = teams[i].points;
+		s0[i].wins = teams[i].wins;
+		s0[i].games = teams[i].games;
+	}
+
+	for (i = 0; i < iterations / 2; i++) {
+		size_t j;
+
+		memcpy(standings, s0, sizeof(standings));
+		memcpy(antistandings, s0, sizeof(antistandings));
+
+		for (j = 0; j < n_games; j++) {
+			size_t wi;
+			size_t li;
+			size_t awi;
+			size_t ali;
+
+			bool antihw;
+			bool antiot;
+			bool hw = urand(games[j].expected, &antihw, seedp);
+			bool ot = urand(otprob, &antiot, seedp);
+
+			wi = games[j].t[!hw];
+			li = games[j].t[hw];
+
+			awi = games[j].t[!antihw];
+			ali = games[j].t[antihw];
+
+			standings[wi].points += 3 - ot;
+			standings[li].points += ot;
+			standings[wi].wins += !ot;
+			standings[wi].games++;
+			standings[li].games++;
+
+			antistandings[awi].points += 3 - ot;
+			antistandings[ali].points += ot;
+			antistandings[awi].wins += !ot;
+			antistandings[awi].games++;
+			antistandings[ali].games++;
+		}
+
+		for (j = 0; j < n_teams; j++) {
+			standings[j].random = getrand(seedp) % (INT_MAX / 2);
+			antistandings[j].random = getrand(seedp) % (INT_MAX / 2);
+		}
+		isort(standings, n_teams, sizeof(*standings), pointscmp);
+
+		for (j = 0; j < n_teams; j++) {
+			teams[standings[j].ti].pointsum += standings[j].points;
+			teams[standings[j].ti].gamessum += standings[j].games;
+			teams[standings[j].ti].poscounts[j]++;
+		}
+
+		isort(antistandings, n_teams, sizeof(*antistandings), pointscmp);
+
+		for (j = 0; j < n_teams; j++) {
+			teams[antistandings[j].ti].pointsum += antistandings[j].points;
+			teams[antistandings[j].ti].gamessum += antistandings[j].games;
+			teams[antistandings[j].ti].poscounts[j]++;
+		}
+	}
+
+}
+
+void simulate_winall(struct team *teams,
+		     size_t n_teams,
+		     const struct game *games,
+		     size_t n_games,
+		     uint32_t otprob,
+		     size_t iterations,
+		     size_t team_id,
+		     bool win,
+		     struct rand_t *seedp)
+{
+	struct team teams_cpy[n_teams];
+	struct game games_cpy[n_games];
+
+	size_t r;
+	size_t w;
+
+	memcpy(teams_cpy, teams, sizeof(teams_cpy));
+
+	for (w = r = 0; r < n_games; r++) {
+		if (games[r].t[0] == team_id || games[r].t[1] == team_id) {
+			size_t opp = games[r].t[0] + games[r].t[1] - team_id;
+			teams_cpy[team_id].points += win * 3;
+			teams_cpy[team_id].games++;
+			teams_cpy[team_id].wins += win;
+
+			teams_cpy[opp].points += !win * 3;
+			teams_cpy[opp].games++;
+			teams_cpy[opp].wins += !win;
+		} else {
+			games_cpy[w++] = games[r];
+		}
+	}
+
+	simulate(teams_cpy, n_teams, games_cpy, w, otprob, iterations, seedp);
+}
+
+void simulate_winloseall(struct team *teams,
+			 size_t n_teams,
+			 size_t team,
+			 const struct game *games,
+			 size_t n_games,
+			 uint32_t otprob,
+			 size_t iterations,
+			 struct rand_t *seedp)
+{
+	struct team teams_cpy[n_teams];
+	size_t poscounts[n_teams][n_teams];
+	size_t e;
+	size_t i;
+
+	memcpy(teams_cpy, teams, sizeof(teams_cpy));
+	memset(poscounts, 0, sizeof(poscounts));
+
+	for (i = 0; i < n_teams; i++)
+		teams_cpy[i].poscounts = poscounts[i];
+
+	simulate_winall(teams_cpy, n_teams,
+			games, n_games,
+			otprob,
+			iterations,
+			team, true,
+			seedp);
+	simulate_winall(teams_cpy, n_teams,
+			games, n_games,
+			otprob,
+			iterations,
+			team, false,
+			seedp);
+
+	for (i = 0; i < n_teams; i++)
+		if (teams_cpy[team].poscounts[i])
+			break;
+	for (e = n_teams; e > 0; e--)
+		if (teams_cpy[team].poscounts[e - 1])
+			break;
+	for (; i < e; i++)
+		teams[team].poscounts[i] = 1;
+}
+
+struct thread_data {
+	struct team *teams;
+	bool *winloseall;
+	size_t n_teams;
+
+	struct game *games;
+	size_t n_games;
+
+	uint32_t otprob;
+
+	size_t iterations;
+	size_t iterated;
+	pthread_mutex_t mutex;
+};
+
+void *simulate_thread(void *data)
+{
+	struct thread_data *td = data;
+	struct team teams_cpy[td->n_teams];
+	size_t i;
+	size_t iters;
+	size_t poscounts[td->n_teams][td->n_teams];
+	size_t n = 0;
+	struct rand_t seedp;
+
+	getseed(&seedp);
+
+	pthread_mutex_lock(&td->mutex);
+	memcpy(teams_cpy, td->teams, sizeof(teams_cpy));
+	pthread_mutex_unlock(&td->mutex);
+
+	memset(poscounts, 0, sizeof(poscounts));
+
+	for (i = 0; i < td->n_teams; i++) {
+		teams_cpy[i].poscounts = poscounts[i];
+		pthread_mutex_lock(&td->mutex);
+		if (!td->winloseall[i]) {
+			td->winloseall[i] = true;
+			pthread_mutex_unlock(&td->mutex);
+
+			simulate_winloseall(teams_cpy, td->n_teams, i,
+					    td->games, td->n_games,
+					    td->otprob,
+					    td->iterations / td->n_teams / 2,
+					    &seedp);
+		} else {
+			pthread_mutex_unlock(&td->mutex);
+		}
+	}
+
+	for (;;) {
+		pthread_mutex_lock(&td->mutex);
+
+		if (td->iterations - td->iterated > 100000)
+			iters = 100000;
+		else
+			iters = td->iterations - td->iterated;
+
+		td->iterated += iters;
+		pthread_mutex_unlock(&td->mutex);
+
+		if (!iters)
+			break;
+		n += iters;
+
+		simulate(teams_cpy, td->n_teams,
+			 td->games, td->n_games,
+			 td->otprob,
+			 iters, &seedp);
+	}
+
+	pthread_mutex_lock(&td->mutex);
+	for (i = 0; i < td->n_teams; i++) {
+		size_t j;
+		td->teams[i].pointsum += teams_cpy[i].pointsum;
+		td->teams[i].gamessum += teams_cpy[i].gamessum;
+
+		for (j = 0; j < td->n_teams; j++)
+			td->teams[i].poscounts[j] += teams_cpy[i].poscounts[j];
+	}
+	pthread_mutex_unlock(&td->mutex);
+
+	return (void *)n;
+}
+
+const char *ltrim(const char *s)
+{
+	while (*s++ == ' ');
+	return s - 1;
+}
+
+int main(int argc, char **argv)
+{
+	int opt;
+	double homeadv = 58;
+	double elok = 32;
+	uint32_t otprob = 0.23 * R_MAX;
+	size_t n;
+	size_t m;
+	size_t i;
+	size_t j;
+	size_t iterations = 1000000;
+	size_t threads = 1;
+
+	while ((opt = getopt(argc, argv, "i:t:a:o:")) != -1) {
+		switch (opt) {
+		case 'i':
+			iterations = strtoull(optarg, NULL, 10);
+			break;
+
+		case 't':
+			threads = strtoull(optarg, NULL, 10);
+			break;
+
+		case 'a':
+			homeadv = strtod(optarg, NULL);
+			break;
+
+		case 'o':
+			otprob = strtod(optarg, NULL) * R_MAX;
+			break;
+
+		case 'k':
+			elok = strtod(optarg, NULL);
+			break;
+
+		default:
+			exit(1);
+			break;
+		}
+	}
+
+	iterations += iterations & 1;
+
+	if (scanf(" %zu", &n) < 1)
+		return 1;
+
+	struct team teams[n];
+	size_t teammap[n];
+	size_t poscounts[n][n];
+	memset(poscounts, 0, sizeof(poscounts));
+
+	for (i = 0; i < n; i++) {
+		if (scanf(" %127s", teams[i].name) < 1)
+			return 1;
+		teams[i].id = i;
+		teams[i].points = 0;
+		teams[i].games = 0;
+		teams[i].wins = 0;
+		teams[i].pointsum = 0;
+		teams[i].gamessum = 0;
+		teams[i].poscounts = poscounts[i];
+		teams[i].elo = 2000;
+	}
+
+	if (scanf(" %zu", &m) < 1)
+		return 1;
+
+	for (i = 0; i < m; i++) {
+		size_t hi, ai;
+		unsigned hs, as;
+		bool ot = false;
+		char det[16] = "";
+
+		switch (scanf(" %zu %zu %u %u%15[^\n]", &hi, &ai, &hs, &as, det)) {
+		case 5:
+			ot = !strcmp(ltrim(det), "OT") ||
+			     !strcmp(ltrim(det), "SO");
+			break;
+
+		case 4:
+			break;
+
+		default:
+			exit(1);
+		}
+
+		double expected = 1. / (1 + pow(10, (teams[ai].elo - teams[hi].elo - homeadv) / 400.0));
+		double actual;
+
+		teams[hi].games++;
+		teams[ai].games++;
+
+		if (hs > as) {
+			if (!ot) {
+				actual = 1;
+				teams[hi].wins++;
+				teams[hi].points += 3;
+			} else {
+				actual = 2./3;
+				teams[hi].points += 2;
+				teams[ai].points += 1;
+			}
+		} else {
+			if (ot) {
+				actual = 1./3;
+				teams[hi].points += 1;
+				teams[ai].points += 2;
+			} else {
+				teams[ai].points += 3;
+				actual = 0;
+			}
+		}
+
+		teams[hi].elo += elok * (actual - expected);
+		teams[ai].elo -= elok * (actual - expected);
+	}
+
+	isort(teams, 1[&teams] - teams, sizeof(*teams), teampointscmp);
+
+	for (i = 0; i < n; i++)
+		teammap[teams[i].id] = i;
+
+	if (scanf(" %zu", &m) < 1)
+		return 1;
+
+	struct game games[m];
+
+	for (i = 0; i < m; i++) {
+		size_t h, a;
+		if (scanf(" %zu %zu", &h, &a) < 2)
+			return 1;
+		games[i].t[0] = teammap[h];
+		games[i].t[1] = teammap[a];
+		games[i].expected = R_MAX / (1 + pow(10.0, (teams[games[i].t[1]].elo  - (teams[games[i].t[0]].elo + homeadv)) / 400.0));
+	}
+
+	if (threads > 1) {
+		pthread_t pthrds[threads];
+		bool winloseall[n];
+		memset(winloseall, 0, sizeof(winloseall));
+		struct thread_data td = {
+			teams, winloseall, n,
+			games, m,
+			otprob,
+			iterations, 0,
+			PTHREAD_MUTEX_INITIALIZER
+		};
+		for (i = 0; i < threads; i++)
+			pthread_create(&pthrds[i], NULL, simulate_thread, &td);
+		iterations = 0;
+		for (i = 0; i < threads; i++) {
+			void *a;
+			pthread_join(pthrds[i], &a);
+			iterations += (size_t)a;
+		}
+	} else {
+		struct rand_t seedp;
+		getseed(&seedp);
+
+		for (i = 0; i < n; i++) {
+			simulate_winloseall(teams, n, i,
+					    games, m,
+					    otprob,
+					    iterations / n / 2,
+					    &seedp);
+		}
+
+		simulate(teams, n,
+			 games, m,
+			 otprob,
+			 iterations,
+			 &seedp);
+	}
+
+	isort(teams, n, sizeof(*teams), teamcmp);
+	for (i = 0; i < n; i++) {
+		printf("%s %u %u %llu %.1f %d", teams[i].name,
+		       teams[i].games, teams[i].points,
+		       teams[i].gamessum / iterations,
+		       (double)teams[i].pointsum / iterations,
+		       (int)round(teams[i].elo));
+		for (j = 0; j < n; j++) {
+			if (teams[i].poscounts[j] == iterations + 1)
+				printf(" +");
+			else if (teams[i].poscounts[j])
+				printf(" %e", (teams[i].poscounts[j] - 1) * 100.0 / iterations);
+			else
+				printf(" -");
+		}
+		printf("\n");
+	}
+
+	return 1;
+}
diff --git a/simuformat b/simuformat
@@ -0,0 +1,92 @@
+#!/usr/bin/perl -w -CSAio
+
+use utf8;
+
+my @seps = @ARGV;
+@ARGV = ();
+
+print <<FOO;
+<html>
+<head>
+<style type="text/css">
+table {
+	border-spacing: 0;
+}
+small {
+	font-size: 50%;
+}
+thead td, thead th {
+	border-bottom: 4px solid black;
+}
+tbody tr > td:nth-child(3) {
+	border-right: 4px solid black;
+}
+td {
+	text-align: center;
+}
+FOO
+if (@seps) {
+	print join(", ", map { "tr:nth-child($_) > *" } @seps);
+	print <<FOO;
+{
+	border-bottom: 1px solid black;
+}
+FOO
+	print join(", ", map { "tr:nth-child(" . ($_ + 1) . ") > *" } @seps);
+	print <<FOO;
+{
+	border-top: 1px solid black;
+}
+FOO
+	print join(", ", map { "tr > *:nth-child(" . ($_ + 6) . ")" } @seps);
+	print <<FOO;
+{
+	border-right: 1px solid black;
+}
+FOO
+	print join(", ", map { "tr > *:nth-child(" . ($_ + 7) . ")" } @seps);
+	print <<FOO;
+{
+	border-left: 1px solid black;
+}
+FOO
+}
+print <<FOO;
+td.no {
+	background: red;
+}
+td.yes {
+	background: cyan;
+	border: 1px solid darkcyan;
+	font-weight: bold;
+}
+td.points {
+	background: inherit;
+}
+td.sep {
+	background: black;
+	width: 3px;
+	height: 3px;
+}
+</style>
+</head>
+<body>
+<table>
+FOO
+my $headed = 0;
+my $j = 0;
+while (<<>>) {
+	chomp;
+	my ($team, $games, $actualpoints, $fingames, $points, $elo, @probs) = split / /;
+	print "<thead><tr><th>Joukkue</th><th>O</th><th>P</th><th>Elo</th><th>O</th><th>Ennuste</th>" . join('', map { '<th>' . $_ . '.</th>' } 1..($#probs + 1)) . "</tr></thead>\n<tbody>\n" unless $headed++;
+	print "<tr><th>$team</th><td class=\"games\">$games</td><td class=\"actualpoints\">$actualpoints</td><td class=\"elo\">$elo</td><td class=\"fingames\">$fingames</td><td class=\"points\">$points</td>" . join("", map { $_ eq '+' ? '<td class="yes">' . $_ . '</td>' : $_ eq '-' ? '<td class="no">' . $_ . '</td>' : sprintf('<td style="border: 1px solid rgb(%.1f%%, 80%%, %.1f%%); padding: -1px; background-color: rgb(%.1f%%, 100%%, %.1f%%)">%s</td>', 80 - 80 * sqrt($_ / 100), 80 - 80 * sqrt($_ / 100), 100 - 100 * sqrt($_ / 100), 100 - 100 * sqrt($_ / 100), ($_ < 0.005 ? $_ > 0 ? sprintf("<small>10<sup>%.0f</sup></small>", log($_) / log(10)) : '<small>ε</small>' : $_ < 9.995 ? sprintf("%.2f", $_) : sprintf("%.1f", $_)) ) } @probs) . "</tr>\n";
+}
+
+my $now = `LC_CTIME=fi_FI.UTF-8 date`;
+print <<FOO;
+</tbody>
+</table>
+<small>Viimeksi päivitetty $now</small>
+</body>
+</html>
+FOO