Glicko 2 Groovy Implementation

Glicko2.groovy

static class Glicko2 {

    double _MU
    double _SIGMA
    double _PHI
    double _TAU
    double _EPSILON


    Glicko2(Double mu, Double sigma, Double phi, Double tau, Double epsilon) {
        this._MU = mu
        this._SIGMA = sigma
        this._PHI = phi
        this._TAU = tau
        this._EPSILON = epsilon
    }

    Glicko2() {
        this._MU = MU
        this._PHI = PHI
        this._SIGMA = SIGMA
        this._TAU = TAU
        this._EPSILON = EPSILON
    }

    Rating create_rating(Double _mu, Double _phi, Double _sigma) {
        double mu = (_mu == null) ? this._MU : _mu
        double phi = (_phi == null) ? this._PHI : _phi
        double sigma = (_sigma == null) ? this._SIGMA : _sigma
        return new Rating(mu, phi, sigma)
    }

    Rating scale_down(Rating rating, Double ratio) {
        ratio = ratio ?: 173.7178
        double mu = ( rating.mu - this._MU ) / ratio
        double phi = rating.phi / ratio
        return create_rating(mu, phi, rating.sigma)
    }

    Rating scale_up(Rating rating, Double ratio) {
        ratio = ratio ?: 173.7178
        double mu = rating.mu * ratio + this._MU
        double phi = rating.phi * ratio
        return create_rating (mu, phi, rating.sigma)
    }

    static Double reduceImpact(Rating rating) {
        return 1.0 / Math.sqrt((1.0 + (3.0 * Math.pow(rating.phi, 2.0)) / Math.pow(Math.PI, 2.0)))
    }

    static Double expectScore(Rating rating, Rating otherRating, Double impact) {
        return (1.0 / (1.0 + Math.exp(-impact * (rating.mu - otherRating.mu))))
    }

    static double phiStar(double sigma, double phi) {
        return Math.sqrt(phi * phi + sigma * sigma)
    }

    Rating volatilize(Rating r) {
        def nR = scale_down(r, null)
        nR.phi = phiStar(nR.sigma, nR.phi)
        return scale_up(nR, null)
    }

    class F {
        double difference
        double variance
        double phi
        double tau
        double alpha

        double f(double x) {
            double tmp = Math.pow(phi, 2.0) + variance + Math.exp(x)
            double _a = Math.exp(x) * (Math.pow(difference, 2.0) - tmp) / (2 * Math.pow(tmp, 2.0))
            double _b = (x - alpha) / Math.pow(tau, 2.0)
            return _a - _b
        }
    }

    static double sign(x) {
        return x <=> 0.0
    }


    def determineSigma(Rating rating, double difference, double variance) {
        """Determines new volatility."""

        double phi = rating.phi
        double difference_squared = Math.pow(difference, 2.0)
        // 1. Let a = ln(s^2), and define f(x)
        double alpha = Math.log(Math.pow(rating.sigma, 2.0))

        def f = new F(difference: difference, variance: variance, phi: phi, tau: this._TAU, alpha: alpha)

        // 2. Set the initial values of the iterative algorithm.
        double a = alpha
        double b = 0.0 // just to stop it whining about potentially unset values
        if (difference_squared > (Math.pow(phi, 2.0) + variance)) {
            b = Math.log(difference_squared - Math.pow(phi, 2.0) - variance)
        }
        else {
            int k = 1
            while (f.f(alpha - k * Math.abs(this._TAU)) < 0) {
                k += 1
                b = alpha - k * Math.abs(this._TAU)
            }
        }
        // 3. Let fA = f(A) and f(B) = f(B)
        double fa = f.f(a)
        double fb = f.f(b)
        // 4. While |B-A| > e, carry out the following steps.
        // (a) Let C = A + (A - B)fA / (fB-fA), and let fC = f(C).
        // (b) If fCfB < 0, then set A <- B and fA <- fB; otherwise, just set
        //     fA <- fA/2.
        // (c) Set B <- C and fB <- fC.
        // (d) Stop if |B-A| <= e. Repeat the above three steps otherwise.
        double c
        double fc
        double d
        double fd
        for (int i = 100; i > 0; i--) {
            if (Math.abs(b-a) <= this._EPSILON) {
                return Math.exp(a / 2.0)
            }

            c = (a + b) * 0.5
            fc = f.f(c)
            d = c + (c-a)*(sign(fa-fb)*fc)/ Math.sqrt(fc*fc-fa*fb)
            fd = f.f(d)

            if (sign(fd) != sign(fc)) {
                a = c
                b = d
                fa = fc
                fb = fd
            } else if (sign(fd) != sign(fa)) {
                b = d
                fb = fd
            } else {
                a = d
                fa = fd
            }
        }

        throw new Exception("Exceeded iterations: Rating: ${rating} / difference: ${difference} / variance: ${variance}")
    }

    Rating rate(Rating _rating, List<Series> series) {
        // Step 2. For each player, convert the rating and RD's onto the
        //         Glicko-2 scale.
        Rating rating = scale_down(_rating, null)
        // Step 3. Compute the quantity v. This is the estimated variance of the
        //         team's/player's rating based only on game outcomes.
        // Step 4. Compute the quantity difference, the estimated improvement in
        //         rating by comparing the pre-period rating to the performance
        //         rating based only on game outcomes.
        double d_square_inv = 0.0
        double variance_inv = 0.0
        double difference = 0.0
        series.each { Series sera ->
            Rating otherRating = scale_down(sera.otherRating, null)
            double  impact = reduceImpact(otherRating)
            double expected_score = expectScore(rating, otherRating, impact)
            variance_inv += Math.pow(impact, 2.0) * expected_score * (1 - expected_score)
            difference += impact * (sera.score - expected_score)
            d_square_inv += (
                    expected_score * (1 - expected_score) *
                            Math.pow(Q, 2.0) * Math.pow(impact, 2.0))
        }
        difference /= variance_inv
        double variance = 1.0 / variance_inv
        double denom = Math.pow(rating.phi , 2.0) + d_square_inv
        double phi = Math.sqrt(1 / denom)
        // Step 5. Determine the new value, Sigma', ot the sigma. This
        //         computation requires iteration.
        double sigma = this.determineSigma(rating, difference, variance)
        // Step 6. Update the rating deviation to the new pre-rating period
        //         value, Phi*.
        double phi_star = Math.sqrt(Math.pow(rating.phi, 2.0) + Math.pow(sigma, 2.0))
        // Step 7. Update the rating and RD to the new values, Mu' and Phi'.
        phi = 1.0 / Math.sqrt(1.0 / Math.pow(phi_star, 2.0) + 1.0 / variance)
        double mu = rating.mu + Math.pow(phi, 2.0) * (difference / variance)
        // Step 8. Convert ratings and RD's back to original scale.

        return scale_up(create_rating(mu, phi, sigma), null)
    }
}