% Script de Matlab para hallar equilibrios de Nash 
% en el problema de las rutas
% Ricard Torres <rtorres@itam.mx>
%
% Num conductores
ncond = 4 ;
% Num individuos en cada tramo, con asignacion = (N_abd, N_acd, N_abcd)
n_ab = @(asign) asign(1) + asign(3);
n_bd = @(asign) asign(1);
n_ac = @(asign) asign(2);
n_bc = @(asign) asign(3); 
n_cd = @(asign) asign(2) + asign(3);
% Tiempos en los distintos tramos
t_ab = @(asign) 2 * n_ab(asign);
t_ac = @(asign) 7 + n_ac(asign);
t_bd = @(asign) 7 + n_bd(asign);
t_bc = @(asign) n_bc(asign);
t_cd = @(asign) 2 * n_cd(asign);
% Tiempos individuales en las distintas rutas (si no son usadas ponemos 0)
t_abd  = @(asign) (t_ab(asign) + t_bd(asign)) * (asign(1) > 0);
t_acd  = @(asign) (t_ac(asign) + t_cd(asign)) * (asign(2) > 0);
t_abcd = @(asign) (t_ab(asign) + t_bc(asign) + t_cd(asign)) * ...
                       (asign(3) > 0);
% Tiempo agregado que corresponde a la asignacion
tagg = @(asign) asign(1) * t_abd(asign) + asign(2) * t_acd(asign) + ...
         asign(3) * t_abcd(asign);
% Vector de tiempos individuales de las diferentes rutas, 
% en formato (T_abd, T_acd, T_abcd)
Vtiemp = @(asign) [ t_abd(asign) , t_acd(asign), t_abcd(asign)];
% Num distintas asignaciones es
%  1 + 2 + ... + (ncond+1) (ya que comenzamos con 0)
nas = (ncond + 1) * (ncond + 2) / 2;
% Listamos en una matriz las distintas asignaciones, 
% en la forma (N_abd, N_acd, N_abcd)
Masign = zeros(nas, 3);
asig = 0;
for i = 0 : ncond
    for j = 0 : ncond-i
      asig = asig + 1;
        Masign(asig,:) = [i, j, ncond - (i+j)] ;
    end
end
% Listamos en una matriz los tiempos individuales
% y en un vector el tiempo agregado para cada asignacion
Mtiemp = zeros(nas, 3);
Vtagg = zeros(nas, 1);
for asig = 1 : nas
     Mtiemp(asig,:) = Vtiemp(Masign(asig,:)) ;
     Vtagg(asig) = tagg(Masign(asig,:)) ;
end
% Determinacion de si una asignacion _no_ es de equilibrio:
% la funcion sera positiva si alguien tiene incentivos a cambiar 
% de ruta, y cero si nadie los tiene
no_equil = @(asign) ...
   (asign(1) > 0) * ...
   (( t_abd(asign) > t_acd([asign(1)-1,asign(2)+1,asign(3)]) ) ...
    + ( t_abd(asign) > t_abcd([asign(1)-1,asign(2),asign(3)+1]) ) ) ...
 + (asign(2) > 0) * ...
   (( t_acd(asign) > t_abd([asign(1)+1,asign(2)-1,asign(3)]) ) ...
    + ( t_acd(asign) > t_abcd([asign(1),asign(2)-1,asign(3)+1]) ) ) ...
 + (asign(3) > 0) * ...
   (( t_abcd(asign) > t_abd([asign(1)+1,asign(2),asign(3)-1]) ) ...
    + ( t_abcd(asign) > t_acd([asign(1),asign(2)+1,asign(3)-1]) ) ) ;
% Tiempo agregado minimo (optimo colectivo):
taggmin = min(Vtagg);
% Miramos que asignaciones son de equilibrio (si hay alguna), y
% cuales son optimas
nash = zeros(nas, 1);
optim = zeros(nas, 1);
for asig = 1 : nas
    nash(asig) = ~ no_equil(Masign(asig,:)) ;
    optim(asig) = ( Vtagg(asig) == taggmin ) ;
end
% Impresion de datos para todas las asignaciones
fprintf('\tN_abd\tN_acd\tN_abcd\tT_abd\tT_acd\tT_abcd\tTAgg\tOptim\tEquil\n');
for asig = 1 : nas
     fprintf('\t%3d\t%3d\t%3d\t%3d\t%3d\t%3d\t%3d', ...
               Masign(asig,:),Mtiemp(asig,:),Vtagg(asig));
     fprintf('\t'); 
     if optim(asig)
        fprintf('%3d',1); 
     end
     fprintf('\t'); 
     if nash(asig)
        fprintf('%3d',1); 
     end
     fprintf('\n');
end
% Impresion de resumen (solo optimos y equilibrios)
fprintf('\nResumen:\n');
fprintf('\tN_abd\tN_acd\tN_abcd\tT_abd\tT_acd\tT_abcd\tTAgg\tOptim\tEquil\n');
for asig = 1 : nas
   if ( nash(asig) + optim(asig) )
     fprintf('\t%3d\t%3d\t%3d\t%3d\t%3d\t%3d\t%3d', ...
         Masign(asig,:),Mtiemp(asig,:),Vtagg(asig));
     if ( nash(asig) + optim(asig) == 2 )
            fprintf('\t%3d\t%3d',1,1); 
     else
         if optim(asig)
            fprintf('\t%3d',1); 
         else
            fprintf('\t\t%3d',1); 
         end
     fprintf('\n');
     end
   end
end

% Output: 
%
%         N_abd   N_acd   N_abcd  T_abd   T_acd   T_abcd  TAgg    Optim   Equil
%           0       0       4       0       0      23      92
%           0       1       3       0      19      20      79
%           0       2       2       0      20      17      74
%           0       3       1       0      21      14      77
%           0       4       0       0      22       0      88
%           1       0       3      19       0      20      79
%           1       1       2      17      17      17      68               1
%           1       2       1      15      18      14      65
%           1       3       0      13      19       0      70
%           2       0       2      20       0      17      74
%           2       1       1      18      15      14      65
%           2       2       0      16      16       0      64       1
%           3       0       1      21       0      14      77
%           3       1       0      19      13       0      70
%           4       0       0      22       0       0      88
% 
% Resumen:
%         N_abd   N_acd   N_abcd  T_abd   T_acd   T_abcd  TAgg    Optim   Equil
%           1       1       2      17      17      17      68               1
%           2       2       0      16      16       0      64       1
%