 < 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

> #start of original code
> norm = require("gaussian")
> norm = norm(0,1)
> 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 parms.pi != undefined
>       @pi = parms["pi"]
>       @tau = parms["tau"]
>     else if parms.mu != undefined
>       @pi = pow(parms["sigma"] , -2)
>       @tau = @pi * parms["mu"]
>     else
>       @pi = 0
>       @tau = 0
>     if isNaN(@pi) || isNaN(@tau) then throw new Error "Gaussian parms can not be NaN"
>   MuSigma:()->
>     ### Return the value of this object as a (mu, sigma) tuple. ###
>     if @pi == 0.0
>       return [0, Infinity]
>     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)->
>     @id = genId()
>     @variables = variables
>     for v in variables
>       v.AttachFactor(this)
>   toString:()->"Factor_#{@id}"
> 
> # 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)->
>     if len(terms_variables) != len(coeffs) then throw new Error "assert error"
>     @sum = sum_variable
>     @terms = terms_variables
>     @coeffs = coeffs
>     super([sum_variable].concat(terms_variables))
>   _InternalUpdate:(variable, y, fy, a)->
>     new_pi = map(range(len(a)), (j)->pow(a[j],2) / (y[j].pi - fy[j].pi))
>     new_pi = 1.0 / sum(new_pi)
>     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}))
>     return
>   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)
>     return
>   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)) when i != index)
>     a.splice(index, 0, 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.apply(this,Array.prototype.slice.call(args, 0))
>     W = @W.apply(this,Array.prototype.slice.call(args, 0))
>     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
> BETA = null
> EPSILON = null
> GAMMA = null
> 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.
>   ###
>   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((a,b)->a.rank-b.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(i.s, new Gaussian({"mu":i.pl.skill[0], "sigma":i.pl.skill[1] + GAMMA})) for i in zip({s:ss, pl:players}))
>   skill_to_perf = (new LikelihoodFactor(i.s, i.p, pow(BETA,2)) for i in zip({s:ss, p:ps}))
>   perf_to_team = (new SumFactor(i.t, [i.p], [1]) for i in zip({p:ps, t:ts}))
>   team_diff = (new SumFactor(i.d, [i.t1, i.t2], [+1, -1]) for i in zip({d:ds, t1:ts.slice(0,-1), t2: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(i.d, (if i.pl1.rank == i.pl2.rank then Vdraw else Vwin), (if i.pl1.rank == i.pl2.rank then Wdraw else Wwin), EPSILON) for i in zip({d:ds, pl1:players.slice(0,-1), pl2: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.
>   #console.log require('util').inspect(ss, true, 10, true)
>   for i in zip({s:ss, pl:players})
>     i.pl.skill = i.s.value.MuSigma()

1,14d0
< #!/usr/bin/python
< # Copyright 2010 Doug Zongker
< #
