% epecgen: generate   s p a r s e   EPEC and write to ampl file.
%
% matlab by Sven Leyffer, Argonne National Laboratory, 2004.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Randomly generate an instance of an EPEC (equilibrium problem
% with equilibrium constraints) and write the resulting model out
% to an AMPL file.
%
% This is an instance of a multi-leader-follower game, with L leaders.
% 
% The problem for leader l=1,...,L is (all have x_l \in R^n):
%
%
%	minimize	x'*G_l*x/2 + g_i'*x
%	subject to	A_l*x_l <= b_l
%		        0 <= y  _|_  N*x + M*y + q >= 0
%
%     where leader l's problem is defined by
%	g_l, G_l objective (quadratic, pos. def.; size is n x L)
%	A_l is p by n random
%	b_l is p by 1 random >= 0 such that A_l*e <= b_l
% 
%     and the followers' problem is defined by
%	N   is m by (n x L) random
%	M   is m by m random \in [0,1] diagonally dominant
%	q   is m by 1 random >= 0
%
%     using a uniform (0,1) random distribution (poss. scaled).
%     To use other distributions replace "rand(" by "randn(" etc.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPECIFY PROBLEM DIMENSIONS, DATA RANGE AND DENSITY HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l          = 4;         % ... number of leaders
p          = 8;	        % ... number of control constraints
m          = 32;        % ... number of states y
n          = 16;        % ... number of controls x
range_A    = [-4, 4];	% ... map [0,1] into [-4,4] here
range_N    = [ 0, 8];	% ... map [0,1] into [ 0,8] here
range_b    = [ 0, 8];	% ... range of values for b
range_g    = [-6, 6];	% ... range of values for g
range_q    = [-4, 4];	% ... range of values for q
dens_G     = 0.02;      % ... density of G_i
dens_A     = 0.2;	% ... density of A
dens_N     = 0.05;	% ... density of N
dens_M     = 0.1;	% ... density of M
outprec    = 0;         % ... output precision (0=2 digits, 1=double)
plotpics   = 1;         % ... provide plots of sparsity pattern or not
%rand('seed',0);         % ... uncomment for different random seed 

% ... precision of output formats for printed outputs
if (outprec) 
   outlenr = '%20.16g  ';  % ... double precision (16 digits)
else
   outlenr = '%10.2g  ';   % ... low precision (2 digits)
end % if
form1i1g = ['\n\t%-6i \t',outlenr];
form2i1g = ['\n\t%-6i %-6i \t',outlenr];
form3i1g = ['\n\t%-6i %-6i %-6i \t',outlenr];

% ... some useful constants
nl = n*l;
e  = ones(n,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATE  S P A R S E   DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ... generate objective Hessian [NB: BLOCK DIAGONAL]
G = zeros(nl,nl,l);
for i=1:l
   % ... generate symmetric random matrix
   G(:,:,i) = sprandsym(nl,dens_G);
   % ... make sure Hessians are positive definite
   mineigG  = eigs(G(:,:,i),1,'SA');
   G(:,:,i) = G(:,:,i) + max(-mineigG,0)*speye(nl,nl);
   % ... find number of nonzeros for plots
   nnzG(i) = nnz( G(:,:,i) );	
end % for
nnzGmax = max( nnzG );           % ... find sparsity pattern of G's
iG      = zeros(nnzGmax,l);
jG      = zeros(nnzGmax,l);
vG      = zeros(nnzGmax,l);
for i=1:l
   [iG(1:nnzG(i),i),jG(1:nnzG(i),i),vG(1:nnzG(i),i)] = find(G(:,:,i));
end % for

% ... generate objective gradient in given range
g = zeros(nl,l);
for i=1:l
   g(:,i) = range_g(1) + rand(nl,1) * (range_g(2) - range_g(1));
end % for

% ... generate leaders' constraint matrices in given range
A = zeros(p,n,l);
b = zeros(p,l);
for i=1:l
   AA         = sprand(p,n,dens_A);
   [iA,jA,vA] = find(AA);
   vA         = range_A(1) + vA * (range_A(2) - range_A(1));
   A(:,:,i)   = full( sparse(iA,jA,vA,p,n) );
   nnzA(i)    = nnz(A(:,:,i));	
   % ... ensure that A*e <= b
   b        = range_b(1) + rand(p,1) * (range_b(2) - range_b(1));
   Ae       = A(:,:,i)*e;
   b(:,i)   = max( b , Ae );
end % for
nnzAmax = max( nnzA );           % ... find sparsity pattern of A's
iA      = zeros(nnzAmax,l);
jA      = zeros(nnzAmax,l);
vA      = zeros(nnzAmax,l);
for i=1:l
   [iA(1:nnzA(i),i),jA(1:nnzA(i),i),vA(1:nnzA(i),i)] = find(A(:,:,i));
end % for

% ... generate follower's LCP: matrix N
N       = sprand(m,nl,dens_N);
[iN,jN] = find(N);
vN      = range_N(1) + N * (range_N(2) - range_N(1));
nnzN    = nnz(N);	
% ... generate follower's LCP: matrix M & ensure it is diagonally dominant
M       = sprand(m,m,dens_M);
dM      = diag(M);				% subtract diagonal from M
MM      = sparse([1:m],[1:m],dM);
Md      = M - MM;			
col_sum = sum(Md,1);				% work out column sum
row_sum = sum(Md,2);				% work out row sum
max_sum = max( col_sum , row_sum' ) + 1E-2;	% max of two above + 1E-2
for i=1:m
   M(i,i) = max_sum(i);				
end % for
nnzM       = nnz(M);	
[iM,jM,vM] = find(M);
% ... generate follower's LCP: right hand side q
q  = range_q(1) + rand(m,1) * (range_q(2) - range_q(1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISPLAY SPARSITY PATTERNS GENERATED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if plotpics
   figure(1); clf; spy(N); title('sparsity pattern of N');
   figure(2); clf; spy(M); title('sparsity pattern of M');
   for i=1:l
      figure(2+2*i-1); clf; spy(A(:,:,i)); title(['sparsity pattern of A ',num2str(i)]);
      figure(2+2*i);   clf; spy(G(:,:,i)); title(['sparsity pattern of G ',num2str(i)]);
   end % for
end % if

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT PROBLEM TO AMPL FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ... open file
fid = fopen('epec.dat','w+');
fprintf (fid,'# random EPEC   s p a r s e   data from matlab\n');

% ... write dimensions of problem
fprintf(fid,'\nparam l := %-4i ;  ',l); fprintf(fid,'# ... number of leaders');
fprintf(fid,'\nparam n := %-4i ;  ',n); fprintf(fid,'# ... number of leaders variables');
fprintf(fid,'\nparam m := %-4i ;  ',m); fprintf(fid,'# ... number of followers');
fprintf(fid,'\nparam p := %-4i ;  ',p); fprintf(fid,'# ... number of leaders constraints');
fprintf(fid,'\n');

% ... Hessian matrices of leaders
fprintf(fid,'\nparam:\t \tG :=');
for i=1:l
   for j=1:nnzG(i)
      fprintf(fid,form3i1g,iG(j,i),jG(j,i),i,vG(j,i));
   end % for
   if (i<l) fprintf(fid,' \n'); end;
end % for
fprintf(fid,' ;\n');

% ... linear objective part of leaders
fprintf(fid,'\nparam:\t \tg :=');
for i=1:l
   for j=1:n
      fprintf(fid,form2i1g,j,i,g(j,i));
   end % for
   if (i<l) fprintf(fid,' \n'); end;
end % for
fprintf(fid,' ;\n');

% ... leaders' constraints matrix A
fprintf(fid,'\nparam:\t \tA :=');
for i=1:l
   for j=1:nnzA(i)
      fprintf(fid,form3i1g,iA(j,i),jA(j,i),i,vA(j,i));
   end % for
   if (i<l) fprintf(fid,' \n'); end;
end % for
fprintf(fid,' ;\n');

% ... leaders' constraints right hand side
fprintf(fid,'\nparam\t: \tb :=');
for i=1:l
   for j=1:p
      fprintf(fid,form2i1g, j, i, b(j,i));
   end % for
end % for
fprintf(fid,' ;\n');

% ... LCP data : N*x 
fprintf(fid,'\nparam:\t \tN :=');
for i=1:nnzN
   fprintf(fid,form2i1g,iN(i),jN(i),vN(i));
end % for
fprintf(fid,' ;\n');

% ... LCP data : M*y 
fprintf(fid,'\nparam:\t \tM :=');
for i=1:nnzM
   fprintf(fid,form2i1g,iM(i),jM(i),vM(i));
end % for
fprintf(fid,' ;\n');

% ... LCP data : q 
fprintf(fid,'\nparam\t: \tq :=');
for j=1:m
   fprintf(fid,form1i1g, j, q(j));
end % for
fprintf(fid,' ;\n');

% ... close AMPL file
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  E N D  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%