 < from scipy.stats.distributions import norm as scipy_norm
 < from math import sqrt
 < norm = scipy_norm()
+< pdf = norm.pdf
+< cdf = norm.cdf
+< icdf = norm.ppf    # inverse CDF
+< # Update rules for approximate marginals for the win and draw cases,
+< # respectively.
+< Vwin=(t, e)->
+<   return pdf(t-e) / cdf(t-e)
+< Wwin=(t, e)->
+<   return Vwin(t, e) * (Vwin(t, e) + t - e)
+< Vdraw=(t, e)->
+<   return (pdf(-e-t) - pdf(e-t)) / (cdf(e-t) - cdf(-e-t))
+< Wdraw=(t, e)->
+<   return pow(Vdraw(t, e) , 2) + ((e-t) * pdf(e-t) + (e+t) * pdf(e+t)) / (cdf(e-t) - cdf(-e-t))
+< class Gaussian
+<   ###
+<   Object representing a gaussian distribution.  Create as:
+<     new Gaussian({mu=..., sigma=...})
+<       or
+<     new Gaussian({pi=..., tau=...})
+<       or
+<     new Gaussian()    # gives 0 mean, infinite sigma
+<   ###
+<   constructor:(parms={})->
+<     if "pi" in parms
+<       @pi = parms["pi"]
+<       @tau = parms["tau"]
+<     else if "mu" in parms
+<       @pi = pow(parms["sigma"] , -2)
+<       @tau = @pi * parms["mu"]
+<     else
+<       @pi = 0
+<       @tau = 0
+<   repr:()->
+<     return "N(pi={0.pi},tau={0.tau})".format(this)
+<   toString:()->
+<     if @pi == 0.0
+<       return "N(mu=0,sigma=inf)"
+<     else
+<       sigma = sqrt(1/@pi)
+<       mu = @tau / @pi
+<       return "N(mu={0:.3f},sigma={1:.3f})".format(mu, sigma)
+<   MuSigma:()->
+<     ### Return the value of this object as a (mu, sigma) tuple. ###
+<     if @pi == 0.0
+<       return 0, float("inf")
+<     else
+<       return @tau / @pi, sqrt(1/@pi)
+<   mul:(other)->
+<     return new Gaussian({"pi":@pi+other.pi, "tau":@tau+other.tau})
+<   div:(other)->
+<     return new Gaussian({"pi":@pi-other.pi, "tau":@tau-other.tau})
+< class Variable
+<   ### A variable node in the factor graph. ###
+<   constructor:()->
+<     @value = new Gaussian()
+<     @factors = {}
+<   AttachFactor:(factor)->
+<     @factors[factor] = new Gaussian()
+<   UpdateMessage:(factor, message)->
+<     old_message = @factors[factor]
+<     @value = @value.div(old_message).mul(message)
+<     @factors[factor] = message
+<   UpdateValue:(factor, value)->
+<     old_message = @factors[factor]
+<     @factors[factor] = value.mul(old_message).div(@value)
+<     @value = value
+<   GetMessage:(factor)->
+<     return @factors[factor]
+< class Factor
+<   ### Base class for a factor node in the factor graph. ###
+<   constructor:(variables)->
+<     @variables = variables
+<     for v in variables
+<       v.AttachFactor(this)
+< # The following Factor classes implement the five update equations
+< # from Table 1 of the Herbrich et al. paper.
+< class PriorFactor extends Factor
+<   ### Connects to a single variable, pushing a fixed (Gaussian) value
+<   to that variable. ###
+<   constructor:(variable, param)->
+<     super([variable])
+<     @param = param
+<   Start:()->
+<     @variables[0].UpdateValue(this, @param)
+< class LikelihoodFactor extends Factor
+<   ### Connects two variables, the value of one being the mean of the
+<   message sent to the other. ###
+<   constructor:(mean_variable, value_variable, variance)->
+<     super([mean_variable, value_variable])
+<     @mean = mean_variable
+<     @value = value_variable
+<     @variance = variance
+<   UpdateValue:()->
+<     ### Update the value after a change in the mean (going "down" in
+<     the TrueSkill factor graph. ###
+<     y = @mean.value
+<     fy = @mean.GetMessage(this)
+<     a = 1.0 / (1.0 + @variance * (y.pi - fy.pi))
+<     @value.UpdateMessage(this, new Gaussian({"pi":a*(y.pi - fy.pi), "tau":a*(y.tau - fy.tau)}))
+<   UpdateMean:()->
+<     ### Update the mean after a change in the value (going "up" in
+<     the TrueSkill factor graph. ###
+<     # Note this is the same as UpdateValue, with @mean and
+<     # @value interchanged.
+<     x = @value.value
+<     fx = @value.GetMessage(this)
+<     a = 1.0 / (1.0 + @variance * (x.pi - fx.pi))
+<     @mean.UpdateMessage(this, new Gaussian({"pi":a*(x.pi - fx.pi),"tau":a*(x.tau - fx.tau)}))
+< class SumFactor extends Factor
+<   ### A factor that connects a sum variable with 1 or more terms,
+<   which are summed after being multiplied by fixed (real)
+<   coefficients. ###
+<   constructor:(sum_variable, terms_variables, coeffs)->
+<     assert len(terms_variables) == len(coeffs)
+<     @sum = sum_variable
+<     @terms = terms_variables
+<     @coeffs = coeffs
+<     super([sum_variable] + terms_variables)
+<   _InternalUpdate:(variable, y, fy, a)->
+<     new_pi = 1.0 / (sum(pow(a[j],2) / (y[j].pi - fy[j].pi) for j in range(len(a))))
+<     new_tau = new_pi * sum(a[j] * (y[j].tau - fy[j].tau) / (y[j].pi - fy[j].pi) for j in range(len(a)))
+<     variable.UpdateMessage(this, new Gaussian({"pi":new_pi, "tau":new_tau}))
+<   UpdateSum:()->
+<     ### Update the sum value ("down" in the factor graph). ###
+<     y = [t.value for t in @terms]
+<     fy = [t.GetMessage(this) for t in @terms]
+<     a = @coeffs
+<     @_InternalUpdate(@sum, y, fy, a)
+<   UpdateTerm:(index)->
+<     ### Update one of the term values ("up" in the factor graph). ###
+<     # Swap the coefficients around to make the term we want to update
+<     # be the 'sum' of the other terms and the factor's sum, eg.,
+<     # change:
+<     #
+<     #    x = y_1 + y_2 + y_3
+<     #
+<     # to
+<     #
+<     #    y_2 = x - y_1 - y_3
+<     #
+<     # then use the same update equation as for UpdateSum.
+<     b = @coeffs
+<     a = [-b[i] / b[index] for i in range(len(b)) if i != index]
+<     a.insert(index, 1.0 / b[index])
+<     v = @terms.slice(0)
+<     v[index] = @sum
+<     y = [i.value for i in v]
+<     fy = [i.GetMessage(this) for i in v]
+<     @_InternalUpdate(@terms[index], y, fy, a)
+< class TruncateFactor extends Factor
+<   ### A factor for (approximately) truncating the team difference
+<   distribution based on a win or a draw (the choice of which is
+<   determined by the functions you pass as V and W). ###
+<   constructor:(variable, V, W, epsilon)->
+<     super([variable])
+<     @var = variable
+<     @V = V
+<     @W = W
+<     @epsilon = epsilon
+<   Update:()->
+<     x = @var.value
+<     fx = @var.GetMessage(this)
+<     c = x.pi - fx.pi
+<     d = x.tau - fx.tau
+<     sqrt_c = sqrt(c)
+<     args = (d / sqrt_c, @epsilon * sqrt_c)
+<     V = @V(*args)
+<     W = @W(*args)
+<     new_val = new Gaussian({"pi":c / (1.0 - W), "tau":(d + sqrt_c * V) / (1.0 - W)})
+<     @var.UpdateValue(this, new_val)
+< DrawProbability=(epsilon, beta, total_players=2)->
+<   ### Compute the draw probability given the draw margin (epsilon). ###
+<   return 2 * cdf(epsilon / (sqrt(total_players) * beta)) - 1
+< DrawMargin=(p, beta, total_players=2)->
+<   ### Compute the draw margin (epsilon) given the draw probability. ###
+<   return icdf((p+1.0)/2) * sqrt(total_players) * beta
+< INITIAL_MU = 25.0
+< INITIAL_SIGMA = INITIAL_MU / 3.0
+< SetParameters=(beta=null, epsilon=null, draw_probability=null, gamma=null)->
+<   ###
+<   Sets three global parameters used in the TrueSkill algorithm.
+<   beta is a measure of how random the game is.  You can think of it as
+<   the difference in skill (mean) needed for the better player to have
+<   an ~80% chance of winning.  A high value means the game is more
+<   random (I need to be *much* better than you to consistently overcome
+<   the randomness of the game and beat you 80% of the time); a low
+<   value is less random (a slight edge in skill is enough to win
+<   consistently).  The default value of beta is half of INITIAL_SIGMA
+<   (the value suggested by the Herbrich et al. paper).
+<   epsilon is a measure of how common draws are.  Instead of specifying
+<   epsilon directly you can pass draw_probability instead (a number
+<   from 0 to 1, saying what fraction of games end in draws), and
+<   epsilon will be determined from that.  The default epsilon
+<   corresponds to a draw probability of 0.1 (10%).  (You should pass a
+<   value for either epsilon or draw_probability, not both.)
+<   gamma is a small amount by which a player's uncertainty (sigma) is
+<   increased prior to the start of each game.  This allows us to
+<   account for skills that vary over time; the effect of old games
+<   on the estimate will slowly disappear unless reinforced by evidence
+<   from new games.
+<   ###
+<   global BETA, EPSILON, GAMMA
+<   if beta is null
+<     BETA = INITIAL_SIGMA / 2.0
+<   else
+<     BETA = beta
+<   if epsilon is null
+<     if draw_probability is null
+<       draw_probability = 0.10
+<     EPSILON = DrawMargin(draw_probability, BETA)
+<   else
+<     EPSILON = epsilon
+<   if gamma is null
+<     GAMMA = INITIAL_SIGMA / 100.0
+<   else
+<     GAMMA = gamma
+< SetParameters()
+< AdjustPlayers=(players)->
+<   ###
+<   Adjust the skills of a list of players.
+<   'players' is a list of player objects, for all the players who
+<   participated in a single game.  A 'player object' is any object with
+<   a "skill" attribute (a (mu, sigma) tuple) and a "rank" attribute.
+<   Lower ranks are better; the lowest rank is the overall winner of the
+<   game.  Equal ranks mean that the two players drew.
+<   This function updates all the "skill" attributes of the player
+<   objects to reflect the outcome of the game.  The input list is not
+<   altered.
+<   ###
+<   players = players.slice(0)
+<   # Sort players by rank, the factor graph will connect adjacent team
+<   # performance variables.
+<   players.sort(key=lambda p: p.rank)
+<   # Create all the variable nodes in the graph.  "Teams" are each a
+<   # single player; there's a one-to-one correspondence between players
+<   # and teams.  (It would be straightforward to make multiplayer
+<   # teams, but it's not needed for my current purposes.)
+<   ss = [new Variable() for p in players]
+<   ps = [new Variable() for p in players]
+<   ts = [new Variable() for p in players]
+<   ds = [new Variable() for p in players.slice(0,-1)]
+<   # Create each layer of factor nodes.  At the top we have priors
+<   # initialized to the player's current skill estimate.
+<   skill = [new PriorFactor(s, new Gaussian({"mu":pl.skill[0], "sigma":pl.skill[1] + GAMMA})) for (s, pl) in zip(ss, players)]
+<   skill_to_perf = [new LikelihoodFactor(s, p, pow(BETA,2)) for (s, p) in zip(ss, ps)]
+<   perf_to_team = [new SumFactor(t, [p], [1]) for (p, t) in zip(ps, ts)]
+<   team_diff = [new SumFactor(d, [t1, t2], [+1, -1]) for (d, t1, t2) in zip(ds, ts.slice(0,-1), ts.slice(1))]
+<   # At the bottom we connect adjacent teams with a 'win' or 'draw'
+<   # factor, as determined by the rank values.
+<   trunc = [new TruncateFactor(d, Vdraw if pl1.rank == pl2.rank else Vwin, Wdraw if pl1.rank == pl2.rank else Wwin, EPSILON) for (d, pl1, pl2) in zip(ds, players.slice(0,-1), players.slice(1))]
+<   # Start evaluating the graph by pushing messages 'down' from the
+<   # priors.
+<   for f in skill
+<     f.Start()
+<   for f in skill_to_perf
+<     f.UpdateValue()
+<   for f in perf_to_team
+<     f.UpdateSum()
+<   # Because the truncation factors are approximate, we iterate,
+<   # adjusting the team performance (t) and team difference (d)
+<   # variables until they converge.  In practice this seems to happen
+<   # very quickly, so I just do a fixed number of iterations.
+<   #
+<   # This order of evaluation is given by the numbered arrows in Figure
+<   # 1 of the Herbrich paper.
+<   for i in range(5)
+<     for f in team_diff
+<       f.UpdateSum()             # arrows (1) and (4)
+<     for f in trunc
+<       f.Update()                # arrows (2) and (5)
+<     for f in team_diff
+<       f.UpdateTerm(0)           # arrows (3) and (6)
+<       f.UpdateTerm(1)
+<   # Now we push messages back up the graph, from the teams back to the
+<   # player skills.
+<   for f in perf_to_team
+<     f.UpdateTerm(0)
+<   for f in skill_to_perf
+<     f.UpdateMean()
+<   # Finally, the players' new skills are the new values of the s
+<   # variables.
+<   for s, pl in zip(ss, players)
+<     pl.skill = s.value.MuSigma()
+< __all__ = ["AdjustPlayers", "SetParameters", "INITIAL_MU", "INITIAL_SIGMA"]
+\ No newline at end of file
 ---
+> ###
+> Implements the player skill estimation algorithm from Herbrich et al.,
+> "TrueSkill(TM): A Bayesian Skill Rating System".
+> ###
 > #define a bunch of python methods for JavaScript
 > sqrt = Math.sqrt
 > pow = Math.pow

< # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
< # See the License for the specific language governing permissions and
< # limitations under the License.
19,27c5,30
< from __future__ import print_function
< __author__ = "Doug Zongker <dougz@isotropic.org>"
< import sys

1,311c1,321
< #!/usr/bin/python
< # Copyright 2010 Doug Zongker
< #
