% genBndQPEC: generate a random bound-constrained QPEC & write to BndQPEC.dat
% matlab function to generate random QPEC instances of the form
%
%  minimize    x^T H x/2 + g^T x
%  subject to  0 <= x1  _|_  x2 >= 0
%              L <= x0 <= U,
%
% where
% x = (x0,x1,x2)     n+2*p dimensional vector of continuous variables
% H                  (n+2*p) x (n+2*p) dimensional sparse Hessian matrix
%                    NB: trailing 2*p x 2*p block of H is positive semi-definite
% g                  n+2*p gradient
% L, U               n finite bounds on x0 variables
%
% User sets parameters for size and range below
%
% Matlab by Sven Leyffer, Argonne National Laboratory, FeBruary, 2020.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% specify problem dimensions and data range here
NumProbs   = 10;           % ... number of distinct instances that are generated
n          = 40;	   % ... number of variables, x_0
p          = 40;  	   % ... number of complementarity constraints
denseH     = 0.25;         % ... density of Hessian
PosSemDef  = 1;            % ... 1=positive semi-definite nxn leading block of H, 0 else
rangeG     = [-10, 10];    % ... range of values for g (uniform distribution)
rangeU     = [ 0,20];	   % ... range of values for U (uniform distribution)
rangeL     = [-10, 10];	   % ... range of values for L (uniform distribution)
% ... choose precision of output below (lower precision is more interesting ;-)
nDigits    = 4;            % ... number of significant digits (adjust output below accordingly)
epsPrec    = 2*10^(-nDigits);
outprec    = ['%',num2str(nDigits+6),'.',num2str(nDigits),'g   '];
outlngt    = ['%',num2str(nDigits+6),'i  '];
% ... uncomment next line to start with same random seed every run
rand('seed',0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ... formats for printed outputs
form1i1g = ['\n\t%-4i \t',outprec];
form1i2g = ['\n\t%-4i \t',outprec,outprec];
form1i5g = ['\n\t%-4i ',outprec,outprec,outprec,outprec,outprec];
form5i   = ['\n\t:',outlngt,outlngt,outlngt,outlngt,outlngt];
form5ie  = ['\n\t:',outlngt,outlngt,outlngt,outlngt,outlngt,'  :='];

% ... full problem dimension
N = n + 2*p;

PSD = '';
if PosSemDef
  PSD='-psd';
else
  PSD='-ind';
end

% ... loop over distinct problem instanecs
for k=1:NumProbs

  % ... generate sparse Hessian matrix H (and make sure it it psd if needed)
  H  = full( sprandsym(N, denseH) );               % ... full matrix for easier printing
  H  = round( H , nDigits);                        % ... round matrix to nDigits
  H1 = H(n+1:N,n+1:N);                             % ... extract trailing 2*p x 2*p block
  MinEig1 = max( 0, -min( eig( H1 ) - epsPrec) );  % ... -epsPrec to ensure H is p.s.d.
  H1 = H1 + round( MinEig1*eye(2*p), nDigits);     % ... make trailing matrix p.s.d.
  H(n+1:N,n+1:N) = H1;
  if (PosSemDef==1)
    MinEig = max( 0, -min( eig( H ) - epsPrec) );  % ... -epsPrec to ensure H is p.s.d.
    H      = H + round( MinEig*eye(N) , nDigits);
  end % if

  % ... generate objective gradient term
  g  = round( rangeG(1) + rand(N,1) * (rangeG(2) - rangeG(1)) , nDigits );

  % ... initialize bounds, and then generate L <= x0 <= U
  L        = zeros(N,1);
  U        = 1E20*ones(N,1);
  U(1:n,1) = round(      rangeU(1) + rand(n,1) * (rangeU(2) - rangeU(1))             , nDigits);
  L(1:n,1) = round( min( rangeL(1) + rand(n,1) * (rangeL(2) - rangeL(1)) , U(1:n,1) ), nDigits);
  
  % ... some useful constants
  NN = floor( N/5 );

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ... open output ampl file
  name = ['BndQPEC-',num2str(n),'-',num2str(p),'-',num2str(denseH),PSD,'-',num2str(k-1),'.dat']
  fid  = fopen(name,'w+');
  fprintf (fid,'# random bound-constrained QPEC data from matlab\n');

  % ... write dimensions of problem
  fprintf(fid,'\nparam n := %-4i ;'  ,n);
  fprintf(fid,'\nparam p := %-4i ;\n',p);

  % ... Hessian matrix H
  fprintf(fid,'\nparam\t \tH');
  for i=1:NN
    ii = (i-1)*5;
    fprintf(fid,form5ie, [ii+1:ii+5]);
    for j=1:N
      fprintf(fid,form1i5g,j,H(j,ii+1:ii+5));
    end % for
    fprintf(fid,'\n');
  end % for
  if ((N - 5*NN) > 0) 
    fprintf(fid,form5i, [5*NN+1:N]); fprintf(fid,' := ');
    for j=1:N
      fprintf(fid,form1i5g,j,H(j,5*NN+1:N));
    end % for
  end % if
  fprintf(fid,' ;\n');
  fprintf(fid,'\nparam\t: \tL, U :=');
  for j=1:n
    fprintf(fid,form1i2g, j, L(j), U(j));
  end % for
  fprintf(fid,' ;\n');
  fprintf(fid,'\nparam\t: \tg :=');
  for j=1:N
    fprintf(fid,form1i1g, j, g(j));
  end % for
  fprintf(fid,' ;\n');
  fclose(fid);
end % for
