Glicko-2 implementation in OCaml

glicko2.ml

(* dump of the unused code written in dec 2022, not thoroughly tested *)

open! Base

(* Implementation of the Glicko-2 algorithm: http://www.glicko.net/glicko/glicko2.pdf

   Written with inspirations from
   https://github.com/lichess-org/lila/blob/668814b/modules/rating/src/main/Glicko.scala
   and https://github.com/TaylorP/Glicko2/blob/c117dea05/glicko/rating.cpp. *)

(* Rating  r, mu
   Rating deviation  RD, phi
   Rating volatility  vol, sigma
   The system constant  tau
   The convergence constant  epsilon *)

type t = { r : float; rd : float; vol : float }

let show { r; rd; vol } =
  Printf.sprintf "R=%f;RD=%f;VOL=%f" r rd vol

type outcome = Loss | Draw | Win

let score_of_outcome = function
  | Loss -> 0.
  | Draw -> 0.5
  | Win -> 1.

type gameresult = outcome * t
(** [(outcome, opponent)] *)

(* Defaults for an unrated player.
   The Glicko-2 paper recommends r=1500, RD=350, vol=0.06.
   Lichess uses r=1500, RD=500, vol=0.09 at the moment. *)
let default = { r = 1500.; rd = 500.; vol = 0.09 }

let min_rating = 600.
let max_rating = 4000.

(* The system constant. The recommended value is between 0.3 and 1.2; Lichess
   uses 0.75. *)
let tau = 0.75

(* The convergence tolerance; 0.000001 is recommended in the paper. *)
let epsilon = 0.000001

(* Rating deviation decrease per day; 0.21436 is used in Lichess *)
let rd_decay_per_day = 0.21436

let max_volatility = 0.1

(* Convergence iteration limit is added as per Lichess' implementation *)
let iteration_limit = 2000

let glicko2_constant = 173.7178

let mu_of_r r = (r -. 1500.) /. glicko2_constant [@@inline]
let phi_of_rd rd = rd /. glicko2_constant [@@inline]
let r_of_mu mu = glicko2_constant *. mu +. 1500. [@@inline]
let rd_of_phi phi = glicko2_constant *. phi [@@inline]

let g phi =
  (* g(phi) = 1 / sqrt(1 + 3 * phi^2 / pi^2) *)
  Float.(let a = phi / pi in 1. / sqrt (1. + 3. * a * a))

let e mu mu' phi' =
  (* E(mu, mu_j, phi_j) = 1 / (1 + exp(-g(phi_j)(mu - mu_j))) *)
  let open Float in
  let exponent = -1. * g phi' * (mu - mu') in
  1. / (1. + exp exponent)

let update { r; rd; vol } (results : gameresult list) =
  assert (not @@ List.is_empty results);
  let open Float in
  (* Convert values to the Glicko-2 scale *)
  let mu = mu_of_r r and phi = phi_of_rd rd and sigma = vol in
  (* v = [\sum_{j=1}^{m} g(phi_j)^2 * E(mu, mu_j, phi_j) * {1 - E(mu, mu_j, phi_j)}]^-1
     delta = v * \sum_{j=1}^{m} g(phi_j) * {s_j - E(mu, mu_j, phi_j)}
     for m opponents *)
  let rec calc v idelta = function
    | (outcome, { r = r'; rd = rd'; _ }) :: xs ->
      let mu' = mu_of_r r' and phi' = phi_of_rd rd' in
      let g = g phi' in
      let e = e mu mu' phi' in
      let s = score_of_outcome outcome in
      let v = v + g * g * e * (1. - e) in
      let idelta = idelta + g * (s - e) in
      calc v idelta xs
    | [] -> v, idelta
  in
  let inv_v, idelta = calc 0. 0. results in
  let v = 1. / inv_v in
  let delta = v * idelta in
  let delta2 = delta ** 2. and phi2 = phi ** 2. in
  let a = log (sigma ** 2.) in
  let f x =
    (*        e^x (delta^2 - phi^2 - v - e^x)     (x - a)
      f(x) = -------------------------------  -  -------
                  2(phi^2 + v + e^x)^2            tau^2 *)
    let exp_x = exp x in
    let num = exp_x * (delta2 - phi2 - v - exp_x) in
    let den = 2. * (phi2 + v + exp_x) ** 2. in
    (num / den) - ((x - a) / tau ** 2.)
  in
  (* The iterative convergence alogrithm *)
  let init_aa = a in
  let init_bb =
    if delta2 > phi2 + v then
      log (delta2 - phi2 - v)
    else
      let rec go k = if f (a - k * tau) < 0. then go (k + 1.) else k in
      a - go 1. * tau
  in
  let rec convergence aa f_a bb f_b i =
    if Int.(i > iteration_limit) then
      failwith "Convergence failure";
    if abs (bb - aa) > epsilon then
      let cc = aa + (aa - bb) * f_a / (f_b - f_a) in
      let f_c = f cc in
      if f_c * f_b < 0. then
        convergence bb f_b cc f_c Int.(i + 1)
      else
        convergence aa (f_a / 2.) cc f_c Int.(i + 1)
    else
      exp (aa / 2.)
  in
  (* Update the volatility *)
  let sigma' = convergence init_aa (f init_aa) init_bb (f init_bb) 0 in
  (* Update the rating and the rating deviation *)
  let phi_star = sqrt (phi2 + sigma' ** 2.) in
  let phi' = 1. / sqrt (1. / phi_star ** 2. + inv_v) in
  let mu' = mu + phi' ** 2. * idelta in
  (* Convert values back to the original scale *)
  let r' = r_of_mu mu' and rd' = rd_of_phi phi' and vol' = sigma' in
  { r = r'; rd = rd'; vol = vol' }

(* TODO: decay by time *)
(* let decay { r; rd; vol } =
  (* phi' = phi_star = sqrt(phi^2 + sigma^2) *)
  let phi = phi_of_rd rd and sigma = vol in
  let phi' = Float.(sqrt (phi ** 2. + sigma ** 2.)) in
  let rd' = rd_of_phi phi' in
  { r; rd = rd'; vol } *)

let%expect_test "example from the paper" =
  (* In the paper, tau=0.5 is used, so we can get slightly different values here *)
  let player = { r = 1500.; rd = 200.; vol = 0.06 } in
  let opponent1 = { r = 1400.; rd = 30.; vol = 0.06 } in
  let opponent2 = { r = 1550.; rd = 100.; vol = 0.06 } in
  let opponent3 = { r = 1700.; rd = 300.; vol = 0.06 } in
  let games = [Win, opponent1; Loss, opponent2; Loss, opponent3] in
  Caml.print_endline @@ show @@ update player games;
  [%expect {| R=1464.050680;RD=151.516504;VOL=0.059991 |}]