% moqpgen: generate   s p a r s e   MOQP and write to ampl file.
%
% matlab by Sven Leyffer, Argonne National Laboratory, 2005.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Randomly generate an instance of an MOQP (multi-objective QP)
% and write the resulting model out to an AMPL file.
%
% The problem has p objective functions, n variables, and m constraints:
%
%	minimize	{ x'*G_k*x/2 + g_k'*x , for k=1,...,p }
%	subject to	A*x >= b
%		          x >= 0  (not generated)
%
%     where 
%	g_k, G_k objective (quadratic, pos. def.; size is n x n)
%	A is m by n random
%	b is m by 1 random >= 0 such that A*e >= b
% 
%     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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p          = 3;	        % ... number of objectives
m          = 20;        % ... number of constraints
n          = 40;        % ... number of variables
range_A    = [-4, 4];	% ... map [0,1] into [-4,4] here
range_b    = [ 0, 8];	% ... range of values for b
range_g    = [-6, 6];	% ... range of values for g
dens_G     = 0.02;      % ... density of G_k
dens_A     = 0.2;	% ... density of A
outprec    = 0;         % ... output precision (0=2 digits, 1=double)
plotpics   = 0;         % ... 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
formatG = ['\nlet G[ %-6i , %-6i , %-6i ] := \t',outlenr,';'];
formatg = ['\nlet g[ %-6i , %-6i ] := \t',outlenr,';'];
formatA = ['\nlet A[ %-6i , %-6i ] := \t',outlenr,';'];
formatb = ['\nlet b[ %-6i ] := \t',outlenr,';'];

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

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

% ... generate objective Hessian [NB: BLOCK DIAGONAL]
G = zeros(n,n,p);
for i=1:p
   % ... generate symmetric random matrix
   G(:,:,i) = sprandsym(n,dens_G);
   % ... make sure Hessians are positive definite
   mineigG  = eigs(G(:,:,i),1,'SA');
   G(:,:,i) = G(:,:,i) + max(-mineigG,0)*speye(n,n);
   % ... 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,p);
jG      = zeros(nnzGmax,p);
vG      = zeros(nnzGmax,p);
for i=1:p
   [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(p,n);
for i=1:p
   g(i,:) = range_g(1) + rand(n,1) * (range_g(2) - range_g(1));
end % for

% ... generate constraint matrices in given range & ensure that A*e >= b
A          = zeros(m,n);
b          = zeros(m,1);
AA         = sprand(m,n,dens_A);
[iA,jA,vA] = find(AA);
vA         = range_A(1) + vA * (range_A(2) - range_A(1));
A          = full( sparse(iA,jA,vA,m,n) );
nnzA       = nnz(A);	
b          = range_b(1) + rand(m,1) * (range_b(2) - range_b(1));
b          = min( b , A*e );


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

if plotpics
   figure(1); clf; spy(A); title('sparsity pattern of A');
   for i=1:p
      figure(i+1); clf; spy(G(:,:,i)); title(['sparsity pattern of G ',num2str(i)]);
   end % for
end % if


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

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

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

fprintf(fid,'# Hessian matrices G\n');
for i=1:p
   for j=1:nnzG(i)
      fprintf(fid,formatG,i,iG(j,i),jG(j,i),vG(j,i));
   end % for
   fprintf(fid,' \n'); 
end % for

fprintf(fid,'# Linear objectives g\n');
for i=1:p
   for j=1:n
      fprintf(fid,formatg,i,j,g(i,j));
   end % for
   fprintf(fid,' \n'); 
end % for

fprintf(fid,'# Linear constraints matrix A\n');
for j=1:nnzA
   fprintf(fid,formatA,iA(j),jA(j),vA(j));
end % for
fprintf(fid,'\n');

fprintf(fid,'# Linear constraints: rhs b\n');
for j=1:m
   fprintf(fid,formatb, j, b(j));
end % for
fprintf(fid,'\n');

fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  E N D  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%