riix.models.skf

  1"""simplified kalman filter"""
  2import math
  3import numpy as np
  4from scipy.stats import norm
  5from riix.core.base import OnlineRatingSystem
  6from riix.utils.math_utils import base_10_sigmoid
  7from riix.utils.math_utils import v_and_w_win_scalar, v_and_w_draw_scalar, norm_cdf
  8
  9LOG10 = math.log(10.0)
 10LOG10_SQUARED = LOG10**2.0
 11
 12
 13def bradley_terry_prob(z):
 14    return base_10_sigmoid(z)
 15
 16
 17def bradley_terry_scalar_prob_g_h(z, outcome):
 18    prob = base_10_sigmoid(z)
 19    g = LOG10 * (outcome - prob)
 20    h = LOG10_SQUARED * prob * (1.0 - prob)
 21    return prob, g, h
 22
 23
 24def thurstone_mosteller_scalar_prob_g_h(z, outcome):
 25    prob = norm_cdf(z)
 26    if outcome != 0.5:
 27        sign_multiplier = (2 * outcome) - 1  # maps 1 to 1 and 0 to -1
 28        v, w = v_and_w_win_scalar(z * sign_multiplier, 0.0)
 29        v = v * sign_multiplier
 30    else:
 31        v, w = v_and_w_draw_scalar(z, 0.0)
 32    return prob, v, w
 33
 34
 35class VSKF(OnlineRatingSystem):
 36    """vector covariance simplified kalman filter"""
 37
 38    rating_dim = 2
 39
 40    def __init__(
 41        self,
 42        competitors: list,
 43        mu_0: float = 0.0,
 44        v_0: float = 1.0,
 45        beta: float = 1.0,
 46        s: float = 1.0,
 47        epsilon: float = 1e-2,
 48        dtype=np.float64,
 49        model='bt',  # legal values are bt, tm, and dd
 50        update_method='online',
 51    ):
 52        super().__init__(competitors)
 53        self.mus = np.zeros(self.num_competitors, dtype=dtype) + mu_0
 54        self.vs = np.zeros(self.num_competitors, dtype=dtype) + v_0
 55        self.has_played = np.zeros(self.num_competitors, dtype=np.bool_)
 56        self.beta = beta
 57        self.beta2 = beta**2.0
 58        self.s = s
 59        self.s2 = s**2.0
 60        self.epsilon = epsilon
 61        self.prev_time_step = 0
 62
 63        if update_method == 'batched':
 64            self.update = self.batched_update
 65        elif update_method == 'online':
 66            self.update = self.online_update
 67        else:
 68            raise ValueError(f'update_method must be one of online or batched, got {update_method}')
 69
 70        if model == 'bt':
 71            self.predict_func = bradley_terry_prob
 72            self.prob_g_h_func = bradley_terry_scalar_prob_g_h
 73        elif model == 'tm':
 74            self.predict_func = norm.cdf
 75            self.prob_g_h_func = thurstone_mosteller_scalar_prob_g_h
 76
 77    def predict(self, matchups: np.ndarray, set_cache: bool = False, **kwargs):
 78        """generate predictions"""
 79        ratings_1 = self.mus[matchups[:, 0]]
 80        ratings_2 = self.mus[matchups[:, 1]]
 81        rating_diffs = ratings_1 - ratings_2
 82        probs = self.predict_func(rating_diffs / self.s)
 83        return probs
 84
 85    @property
 86    def ratings(self):
 87        return self.mus
 88
 89    def get_pre_match_ratings(self, matchups: np.ndarray, **kwargs):
 90        mus = self.mus[matchups]
 91        vs = self.vs[matchups]
 92        ratings = np.concatenate((mus[..., None], vs[..., None]), axis=2).reshape(mus.shape[0], -1)
 93        return ratings
 94
 95    def time_dynamics_update(self, time_step, matchups):
 96        """called once per period to model the increase in variance over time"""
 97        active_in_period = np.unique(matchups)
 98        self.has_played[active_in_period] = True
 99        time_delta = time_step - self.prev_time_step
100
101        # update parameters for passage of time
102        beta_t = self.beta**time_delta
103        self.mus[self.has_played] = beta_t * self.mus[self.has_played]
104        self.vs[self.has_played] = (beta_t**2) * self.vs[self.has_played] + (time_delta * self.epsilon)
105        self.prev_time_step = time_step
106
107        return active_in_period
108
109    def batched_update(self, matchups, outcomes, time_step, use_cache=False):
110        """apply one update based on all of the results of the rating period"""
111        active_in_period = self.time_dynamics_update(time_step, matchups)
112        masks = np.equal(matchups[:, :, None], active_in_period[None, :])  # N x 2 x active
113
114        mus = self.mus[matchups]
115        vs = self.vs[matchups]
116
117        omegas = np.sum(vs, axis=1)
118        probs = base_10_sigmoid((mus[:, 0] - mus[:, 1]) / self.s)
119
120        g = LOG10 * (outcomes - probs)
121        h = LOG10_SQUARED * probs * (1.0 - probs)
122
123        denom = (self.s2) + h * omegas
124        mu_updates = (self.s * g) / denom
125        v_updates = h / denom
126
127        mu_updates = mu_updates[:, None].repeat(2, 1)
128        mu_updates[:, 1] *= -1
129
130        v_updates = 1.0 - vs * v_updates[:, None]
131
132        mu_updates_pooled = (mu_updates[:, :, None] * masks).sum(axis=(0, 1))
133        v_updates_pooled = (v_updates[:, :, None] * masks).max(axis=(0, 1))
134
135        self.mus[active_in_period] += self.vs[active_in_period] * mu_updates_pooled
136        self.vs[active_in_period] *= v_updates_pooled
137
138    def online_update(self, matchups, outcomes, time_step, **kwargs):
139        """treat the matchups in the rating period as if they were sequential"""
140        self.time_dynamics_update(time_step, matchups)
141        for idx in range(matchups.shape[0]):
142            comp_1, comp_2 = matchups[idx]
143            mu_1 = self.mus[comp_1]
144            mu_2 = self.mus[comp_2]
145            v_1 = self.vs[comp_1]
146            v_2 = self.vs[comp_2]
147            outcome = outcomes[idx]
148
149            omega = v_1 + v_2
150            z = (mu_1 - mu_2) / self.s
151
152            prob, g, h = self.prob_g_h_func(z, outcome)
153
154            denom = (self.s2) + h * omega
155            mu_update = (self.s * g) / denom
156
157            v_update = h / denom
158
159            self.mus[comp_1] += v_1 * mu_update
160            self.mus[comp_2] -= v_2 * mu_update
161            self.vs[comp_1] *= 1.0 - v_1 * v_update
162            self.vs[comp_2] *= 1.0 - v_2 * v_update
163
164    def print_leaderboard(self, num_places):
165        sort_array = self.mus - (3.0 * np.sqrt(self.vs))
166        sorted_idxs = np.argsort(-sort_array)[:num_places]
167        max_len = min(np.max([len(comp) for comp in self.competitors] + [10]), 25)
168        print(f'{"competitor": <{max_len}}\t{"mu - (3 * sd)"}')
169        for p_idx in range(num_places):
170            comp_idx = sorted_idxs[p_idx]
171            print(f'{self.competitors[comp_idx]: <{max_len}}\t{sort_array[comp_idx]:.6f}')