% Created on June 26, 2008 by Paul Tseng
% Last updated October 14, 2008
%
% This Matlab code implements two parametrized VI approaches to 
% generate solutions of linear-quadratic GNEP.
% Reference: K. Nabetani, P. Tseng, and M. Fukushima,
% Parametrized Variational Inequality Approaches to Generalized Nash 
% Equilibrium Problems with Shared Constraints October 2008.
%
% Input:  N, m, l
%         n_1,..., n_N
%         n by n 		Q 	(n=n_1+...+n_N)
%         n by 1        q
%         m by n        B 
%         m by 1        b
%         l by n		A		
%	    l by 1		a
%         n by 1  	u	
%
%with l=l_1+...+l_N, A = [ A_1 0 ... 0 ; 0 A_2 0 .. 0 ; ... ], a = [ a_1 ; a_2 ; ...]
%B = [ B_1  B_2 ... B_N ].
%See Subsection 3.4.3 of the above reference for the meaning of these notations.
%
%IMPORTANT NOTE: The upper bound u is used in determining the sampling range 
%for resource-direc. parametrized VI, so this bound should be as sharp 
%as possible to limit the sampling range.
%
%Each parametrized VI is formulated as a mixed LCP and solved by 
%pathlcp of (Ferris and Munson).  Other LCP solvers can also be used. 
%
%IMPORTANT NOTE: For resource-direc. param. VI, the mixed LCP may have no solution
%and this should be flagged by setting z(1)<0.  If you use pathlcp, then make
%the following changes to the file pathlcp.m near the end:
%
%    if (status ~= 1)
% %%   status,
% %%   error('Path fails to solve problem');
%     z(1)=-1; 
%    end

clear;
prob=input(' enter benchmark problem # (1 or 2 or 3) ');

if prob==1
%Reference: P. T. Harker, Generalized {Nash} games and quasi-variational
%inequalities, European Journal of Operational Research, 54 (1991), pp.~81--94.
  N=2;
  m=1;
  l=0;
  n=[1 1];
  Q=[2 8/3 ; 5/4 2];
  q=[-34 ; -24.25];
  B=[1 1];
  b=15;
  A=[];
  a=[];
  u=[10 ; 10];
  fprintf(' 1. Harker example, N= %g  n= %g  m = %g \n',N,sum(n),m);
end
if prob==2
%Reference: J.~B. Krawczyk and S.~Uryasev}, Relaxation algorithms to find {Nash}
%equilibria with economic applications, Environmental Modeling and
%Assessment, 5 (2000), pp.~63--73.
  N=3;
  m=2;
  l=0;
  n=[1 1 1];
  Q=[.04 .01 .01 ; .01 .12 .01 ; .01 .01 .04];
  q=[-2.9 ; -2.88 ; -2.85];
  B=[3.25 1.25 4.125 ; 2.2915 1.5625 2.8125];
  b=[100 ; 100];
  A=[];
  a=[];
  u=ones(sum(n),1)*100; 	%upper bound deduced from Bx <= b
  fprintf(' 2. River basin pollution game, N= %g  n= %g  m = %g \n',N,sum(n),m);
end
if prob==3
%Reference: J.-S. Pang and M.~Fukushima,  Quasi-variational inequalities,
%generalized {Nash} equilibria, and multi-leader-follower games,
%Computational Management Science, 2 (2005), pp.~21--56 and Erratum.  
%Modification: shipping cost e_{ij}=0 (instead of e_{ij}=1) for all i=j.
%
%x=(x_{1,11}, x_{1,12}, x_{1,13}, x_{1,21}, x_{1,22}, x_{1,23}, x_{2,21}, x_{2,22}, x_{2,23}, x_{2,31}, x_{2,32}, x_{2,33})
%S=[ E1 E2 ]*x,   p=P0-(P0./Q0).*S  
  N=2;
  m=6;
  l=4;
  n=[6 6];
  P0=[ 40 35 32 ]';
  Q0=[ 500 400 600 ]';
  E=[ 1 1 -1 0 -1 0 ; -1 0 1 1 0 -1 ; 0 -1 0 -1 1 1 ];  	%node-arc incid. matrix
  E1=[1 0 0 1 0 0; 0 1 0 0 1 0 ; 0 0 1 0 0 1 ];
  E2=E1;
  D=diag(P0./Q0);
  Q=[ 2*E1'*D*E1  E1'*D*E2 ; E2'*D*E1 2*E2'*D*E2 ];
  q=15*ones(sum(n),1) + [0 1 1 1 0 1 1 0 1 1 1 0]' - [ E1'*P0 ; E2'*P0 ];  
%  q=15*ones(sum(n),1) + [1 1 1 1 1 1 1 1 1 1 1 1]' - [ E1'*P0 ; E2'*P0 ];  %Fukushiam-Pang

  B=E'*D*[ E1  E2 ];
  b=ones(m,1)+E'*P0;
  A=[ 1 1 1 0 0 0  0 0 0 0 0 0 ; ...
      0 0 0 1 1 1  0 0 0 0 0 0 ;  ...
      0 0 0 0 0 0  1 1 1 0 0 0 ;  ...
      0 0 0 0 0 0  0 0 0 1 1 1 ];
  a=[ 100 ; 50 ; 100 ; 50 ];
  u=[ 100 ; 100 ; 100  ; 50 ; 50 ; 50 ; 100 ; 100 ; 100 ; 50 ; 50 ; 50 ];
  fprintf(' 3. Electric power market, N= %g  n= %g  m = %g \n',N,sum(n),m);
end

sumn=sum(n);
nsamples=input(' enter no. samples per interval (e.g., 20) ');
pvichoice=input(' enter 1 for price- and 2 for resource-dir. param. VI ');
samplechoice=input(' enter 1 for random and 2 for grid sampling ');
printchoice=input(' enter 1 to print detailed intermediate results (else 0, say) ');
tol=1e-6; 		%tolerance for solving parametrized VI and checking for GNE.
infty=1e30;		%set a large value for infinity.
maxtrial=200;	%max # trials per subset K without finding GNE (in price-direc. param VI).

if pvichoice==1

%%%%%%% Price-directed parameterized VI approach %%%%%%%%%

rho=input(' enter rho, where [0,rho] is sampling interval (e.g., 2, 20) ');
%larger rho means more time, but searches over greater parameter space
maxKsize=input(' enter max limit on |K| <= m (e.g., 2) ');
%larger maxKsize means more time, but searches over more subsets K

t0=cputime;

M=[ Q  B' A' ; [-B ; -A] zeros(m+l,m+l) ];
d=[ u ; ones(m+l,1)*infty ];  
nM=size(M,1);

ngne=0;
totalsample=0;
ngne_dim1=0;
for vi_num=0:(N+1)^m-1
  num=vi_num;
  sigma=zeros(m,1);
  for i=m-1:-1:0
    if ( floor(num/(N+1)^i) ~= 0 )
      sigma(i+1)=floor(num/(N+1)^i);
      num=num - sigma(i+1)*(N+1)^i;
    end 
  end
  sigma_map=sigma'

  indx=find(sigma~=0);
  indx0=find(sigma==0);
  Ksize=size(indx,1);
  nintervals=(N-1)*Ksize;
  gridsz=rho/nsamples;
  ntrials=nsamples^nintervals;

  if Ksize<=maxKsize

%pause

    fprintf(' |K|= %g  #trials= %g  \n',Ksize,ntrials);
    flag=0;  	%set to 1 if a gne is found for this sigma
    for trial=1:ntrials

      if samplechoice==1
%%%% random sampling %%%%%
        w=rand(m,N)*rho;
        w(indx0,:)=0;
        for i=1:m
          if sigma(i)~=0 
            w(i,sigma(i))=0;
          end
        end
      else
%%%% grid sampling %%%%%
        num=trial-1;
        gridnum=zeros(nintervals,1);
        for i=nintervals-1:-1:0
          if ( floor(num/nsamples^i) ~= 0 )
            gridnum(i+1)=floor(num/nsamples^i);
            num=num - gridnum(i+1)*nsamples^i;
          end 
        end
%gridnum(i) is the grid number (from 0 to nsamples-1) for the ith interval 
%[0,rho], i=1,...,nintervals.
        w=zeros(m,N);
        intervalnum=0;
        for i=1:m
          if sigma(i)~=0 
            for j=1:N
              if j~=sigma(i)
                intervalnum=intervalnum+1;
                w(i,j)=(gridnum(intervalnum)+1)*gridsz;
              end
            end
          end
        end
      end
      totalsample=totalsample+1;

% Compute q in parameterized VI
      qw=q;
      for j=1:N
        n1=sum(n(1:j-1))+1;
        n2=n1+n(j)-1;
        qw(n1:n2)=qw(n1:n2)+B(:,n1:n2)'*w(:,j);
      end
    
% solve VI(Qx+c,{x|Bx<=b,Ax<=a}).  First transform it into a mixed LCP
% 0 = mid{z, M*z+c, z-d} or, equivalently, z = min(max(z-(M*z+c),0),d).
      c=[ qw ; b ; a ];
      z=pathlcp(M,c,zeros(nM,1),d);
      x=z(1:sumn);

% Check if solution of LCP is a GNE
      err=max(abs((B*x-b)'*w));
      if err<tol 
        flag=1;
        if ngne>0 & norm(gne(1:nM,ngne)-z)<1e-10
          if printchoice==1 fprintf(' same GNE found! \n'); end
        else 
          ngne=ngne+1;
          if Ksize==1
            ngne_dim1=ngne_dim1+1;
          end
          gne(1:nM,ngne)=z;
          if printchoice==1 fprintf(' maybe new GNE found! \n'); end
        end
      else
        if printchoice==1 fprintf(' GNE not found,  err= %g \n',err); end
      end
      if flag==0 & trial>maxtrial
        fprintf(' no GNE found for this K after %g trials, exit \n',trial);
        break;
      end
    end

  end
  fprintf(' #GNEs found so far = %g  \n',ngne);

end

fprintf(' #LCPs solved= %g   cpu= %g \n',totalsample,cputime-t0);
fprintf(' #GNEs found from sampling over those K with |K|=1 = %g \n',ngne_dim1);


else
%%%%%%% Resource-directed parameterized VI approach %%%%%%%%%

rho=input(' enter rho>0, where -rho is a lower bd on beta (e.g., 10) ');

t0=cputime;
for i=1:m
  gmin=0;
  for j=1:N
    n1=sum(n(1:j-1))+1;
    n2=n1+n(j)-1; 
    gmin=min(gmin,min(B(i,n1:n2),0)*u(n1:n2));
  end
  t=-b(i)/N + gmin;
  alpha(i)=max(t,-rho);
end

%Sample (beta_{i,1},...,beta_{i,N}) from the simplex with vertices
%(1-N,1,...,1), (1,1-N,...,1), ..., (1,...,1,1-N), and multiply by max(alpha,-rho)

ngne=0;
if samplechoice==1
  if N==2  ntrials=nsamples^m;  end
  if N==3  ntrials=(nsamples*(nsamples+1)/2)^m; end
  if N>3  ntrials=nsamples^((N-1)*m);  end	%replace this by exact formula?
else
  nintervals=(N-1)*m;
  gridsz=1/(nsamples-1);
  ntrials=nsamples^nintervals;
end
fprintf(' #trials= %g \n',ntrials);
trial=0;
totalsample=0;
while (trial<ntrials)

  trial=trial+1;    

  if samplechoice==2
%%%% grid sampling %%%%%
    num=trial-1;
    gridnum=zeros(nintervals,1);
    for i=nintervals-1:-1:0
      if ( floor(num/nsamples^i) ~= 0 )
        gridnum(i+1)=floor(num/nsamples^i);
        num=num - gridnum(i+1)*nsamples^i;
      end 
    end
%gridnum(i) is the grid number (from 0 to nsamples-1) for dimension i, 
%i=1,...,nintervals.
    intervalnum=0;
  end

  acceptrial=1;
  for i=1:m
    if samplechoice==1
%%%% random sampling %%%%%
      wt=rand(1,N);
      wt=wt/sum(wt);
    else
%%%% grid sampling %%%%%
      wt=zeros(1,N);		%needed?
      for j=1:N-1
        intervalnum=intervalnum+1;
        wt(j)=gridnum(intervalnum)*gridsz;
      end
      if sum(wt(1:N-1))<=1
        wt(N)=1-sum(wt(1:N-1));
      else
        acceptrial=0;		%trial corresp. unusable weight
        break;
      end
    end
    beta(i,1:N)=alpha(i)*(ones(1,N)-N*wt);
  end

%pause;

  if acceptrial==1

    totalsample=totalsample+1;
    mm=m*N;
    BB=zeros(mm,sumn);
    bb=zeros(mm,1);
    for j=1:N
      n1=sum(n(1:j-1))+1;
      n2=n1+n(j)-1;
      BB((j-1)*m+1:j*m,n1:n2)=B(:,n1:n2);
      bb((j-1)*m+1:j*m)=beta(:,j)+b/N;
    end

% solve VI(Qx+c,{x|BBx<=bb,Ax<=a}).  First transform it into a mixed LCP
% 0 = mid{z, M*z+c, z-d} or, equivalently, z = min(max(z-(M*z+c),0),d).
    M=[ Q  BB' A' ; [-BB ; -A] zeros(mm+l,mm+l) ];
    c=[ q ; bb ; a ];
    d=[ u ; ones(mm+l,1)*infty ];
    
    nM=size(M,1);
    
    clear z;
% Must first modify pathlcp.m to set z(1)<0 when Path fails to solve LCP.
    z=pathlcp(M,c,zeros(nM,1),d);

    if z(1)>=0 
% Check if solution of LCP is a GNE
      x=z(1:sumn);
      lambda=z(sumn+1:sumn+mm);
      s=BB*x-bb;
      gneflag=1;
      for i=1:m
        indx=[0:N-1]*m+i;
        err1=max(abs(s(indx)));
% Check 
%        err2=max(abs(lambda(indx)));
%        if err1>tol & err2>tol
        err3=max(s(indx));
        if err1>tol & err3>-tol
          gneflag=0;
          break;
        end
      end
      if gneflag==1
        if ngne>0 & norm(gne(1:nM,ngne)-z)<1e-10
          if printchoice==1 fprintf(' same GNE found! \n'); end
        else 
          if printchoice==1 fprintf(' maybe new GNE found! \n'); end
        end
        ngne=ngne+1;
        gne(1:nM,ngne)=z;
      else
        if printchoice==1 fprintf(' GNE not found,  err1= %g  err2= %g \n',err1,err2); end
      end

    else
      if printchoice==1 fprintf(' parametrized VI likely infeasible! \n'); end
    end
    if mod(totalsample,100)==0  
       fprintf(' #GNEs found after %g trials = %g \n',totalsample,ngne);
    end

  end
end

fprintf(' #LCPs solved= %g  cpu= %g \n',totalsample,cputime-t0);

end

%%% Check gnes for duplicates %%%
ngne1=min(ngne,1);
for k=2:ngne
   err=min(max(abs(gne(1:sumn,1:k-1)-gne(1:sumn,k)*ones(1,k-1))));
   if err>tol
     ngne1=ngne1+1;
   end
end
fprintf(' #GNEs found= %g,  #distinct GNEs (1-norm dist. > %g) = %g \n',ngne,tol,ngne1);