%% Ranking teams with linear systems.
% Here, we'll apply an idea from the ranking system:
% http://www.raymondcheong.com/rankings/background.html
% Later in the class, we'll see that this is a special case
% of another ranking system, the Massey method.
% This will also be the first time we see the Laplacian matrix,
% but we'l derive the method without using it.

%% The ranking method
% Copied from http://www.raymondcheong.com/rankings/background.html
% A team's ASM depends on the strength of the opponents. 
% Suppose we had a measure of each team's strength, $\vr$ 
% (i.e. the ranking), such that $r_i - r_j$ is the expected 
% score margin when team $i$ plays team $j$. For $\vr$ to be 
% consistent with the available data, then the average 
% discrepancy between the actual and expected score margin
% should be zero.  In other words, for each team $i$, $\vr$
% ought to satisfy the equation:
%
% \[ \sum_{j \in Games(i)} (r_i - r_j) - (S_i - S_j) = 0 \]
%
% where Games(i) is the set of other teams that team $i$ played,
% and S_i and S_j are the scores of team $i$ and $j$ when they play.

%% Constructing the linear system
%
% In this equation, the data are S_i, the scores.  The model is
% \sum_{j \in Games(i)} (r_i - r_j).  The variables are $r_i$,
% the unknown "ranks".  We can build them into a linear system
% Ax = b
% where b_i = sum_j \in Games(i) (S_i - S_j)
% and   x_i = r_i
% and   \sum_{j} A_{i,j} x_j = \sum_{j \in Games(i)} (r_i - r_j).
%
% We build the matrix A, by iterating over all games and putting
% in a (1,-1) pair in the appropriate row.
%
% [Explain this on the board]

%% A real example
% Here, the data comes from Amy Langville's forthcoming book
% on ranking problems.  

teams = {'duke';'miami';'unc';'uva';'vt'};
data = [ % team 1 team 2, team 1 pts, team 2 pts
    1 2  7 52;
    1 3 21 24;
    1 4  7 38;
    1 5  0 45;
    2 3 34 16;
    2 4 25 17;
    2 5 27  7;
    3 4  7  5;
    3 5  3 30;
    4 5 14 52
];

%% Building the system

% dimension of the system
ngames = size(data,1);
nteams = length(teams);

% construct the empty variables
A = zeros(nteams, nteams);
b = zeros(nteams, 1);

% loop over all games
for g=1:ngames
  i = data(g,1);
  j = data(g,2);
  Si = data(g,3);
  Sj = data(g,4);
  
  % work from team i's perspective
  A(i,i) = A(i,i) + 1;
  A(i,j) = A(i,j) - 1;
  b(i) = b(i) + Si - Sj;

  % now from team j's perspective.  
  A(j,j) = A(j,j) + 1;
  A(j,i) = A(j,i) - 1;
  b(j) = b(i) + Sj - Si;
end

%% Solving the system
% For this class, we can use Matlab's backslash.
%
% [In class discussion of Backslash and inv.]

x = A\b

%%
% That answer makes no sense.  Let's check the residual
r = A*x - b

%%
% Hmm... it isn't zero, what could be the problem?
% Later, we'll see a better way of checking this, but
% let's make sure it's non-singular.
rank(A)

%% 
% Ahah, the rank is 4, but it should have been 5 for a non-singular
% system.  
%
% How could we have known that this system was singular?
%
% [Discussion of original equation.]

%% Fixing the singularity.
% The issue is that we can add a constant to any "rank" r_i.
% So let's "fix" these scores around a particular point, and
% see if that will eliminate the singularity.  What we'll do
% is just set r_1 = 1.  We can do this by changing the system
% of equations to "set" r_1 to 1 as follows.

Aorig = A; % save for later
A(1,:) = 0;
A(1,1) = 1;
borig = b;
b(1) = 1;

%% 
% Did it work?
rank(A)

%%
% What are the ranks now?
x = A\b

%%
% Check the residual
Aorig*x - borig
%% 
% Those look better!  Now, which team was the best?
teams

%% Discussion
% Later on in the class we'll see how this is actually an instance
% of the "Massey" method for ranking and how the matrix A here is
% really a Laplacian matrix for the complete graph:

Aorig

