%% Ranking teams with least squares
% This example shows how to construct a least squares problem to determine
% the strength of a sport's team.  In the particular data we have, it'll be
% a basketball team.

%% The model
% Our model is quite simple.  We pick one team, and call them the reference
% team.  The reference team will have a score of "0".  Whenever
% two teams play a game, we expect the score difference of the game to be the 
% difference in their strengths:
%  points(i) - points(j) = strength(i) - strength(j)
% Note that, when written for all teams, this is a matrix equation:
%  p = G*s
% where
%  p(k) = points(i) - points(j) (i.e. the kth game is between teams i and j)
%  G(k,i) = 1 for team i, -1 for team j (in the kth game) and 0 otherwise.
%    in the case that one of the teams is the reference team, just ignore
%    that entry.
%  s(i) = strength(i)
% Of course, this won't hold exactly, because teams don't "exactly" obey
% this model.  However, if we minimize the sum of squares of the difference
% between the predictions for the games, and what actually happens, we
% arrive at a least-squares problem.


%% 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
G = zeros(ngames, nteams-1);
p = zeros(ngames, 1);
refteam = nteams;

% loop over all games
for g=1:ngames
  i = data(g,1);
  j = data(g,2);
  Si = data(g,3);
  Sj = data(g,4);
  
  if i~=refteam
      G(g,i) = 1;
  end
  if j~=refteam
      G(g,j) = -1;
  end
  p(g) = Si - Sj;
end

%% Solving the system
% For this class, we can use Matlab's backslash, which knows how to solve
% least squares problems

spartial = G\p
sfull = [spartial; 0]

%% How close are we?
predicted = G*spartial;

for g=1:ngames
    fprintf('%10s vs. %10s: predicted = %4.1f   actual = %4.1f\n', ...
        teams{data(g,1)}, teams{data(g,2)}, predicted(g), p(g));
end