% write_qpec: read QPEC from .mat file and write to ampl file
% 
% (1) load a qpec of type 300 (LCP-QP) from qpecgen_data.mat
%
%	min	0.5*(x,y)'*P*(x,y) + c'x + d'y
%	s.t.	A*(x,y) + a <= 0		(1)
%		F(x,y) = N*x + M*y + q >= 0
%		y >= 0
%		F(x,y)*y = 0
%
% (2) write data to qpec.dat file for AMPL
%
% Sven Leyffer, University of Dundee, May 2000.

% ... set output precision here
outprec = '%20.16g  ';                   % double precision (16 digits)
outlngt = '%20i  ';
%outprec = '%12.4g  ';                    % low precision (4 digits)
%outlngt = '%12i  ';

% ... formats for printed outputs
form1i1g = ['\n\t%-4i \t',outprec];
form1i3g = ['\n\t%-4i ',outprec,outprec,outprec];
form3i   = ['\n\t:',outlngt,outlngt,outlngt];
form3ie  = ['\n\t:',outlngt,outlngt,outlngt,'  :='];

% .. read data from file generated by qpecgen
load qpecgen_data;

% .. problem dimensions
n_x = length(c);
n_y = length(d);
m_1 = length(a);

% ... some useful constants
mm = floor( n_y/3 );
nn = floor( n_x/3 );

% ... open output ampl file 
fid = fopen('qpec.dat','w+');
fprintf (fid,'# QPEC data from qpecgen & matlab\n');

% ... write dimensions
fprintf(fid,'\nparam n_x := %-4i ;',n_x);
fprintf(fid,'\nparam n_y := %-4i ;',n_y);
fprintf(fid,'\nparam m_1 := %-4i ;\n',m_1);

% ... linear part of objective c'x & d'y
fprintf(fid,'\nparam\t: \tc :=');
for i=1:n_x
   fprintf(fid,form1i1g, i, c(i));
end % for
fprintf(fid,' ;\n');
fprintf(fid,'\nparam\t: \td :=');
for i=1:n_y
   fprintf(fid,form1i1g, i, d(i));
end % for
fprintf(fid,' ;\n');

% ... quadratic part of objective Pxx
fprintf(fid,'\nparam\t \tPxx');
for i=1:nn
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:n_x
      fprintf(fid,form1i3g,j,P(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_x - 3*nn) > 0) 
   fprintf(fid,form3i, [3*nn+1:n_x]); 
   fprintf(fid,' := ');
   for j=1:n_x
      fprintf(fid,form1i3g,j,P(j,3*nn+1:n_x));
   end % for
end % if
fprintf(fid,' ;\n');

% ... quadratic part of objective Pxy
fprintf(fid,'\nparam\t \tPxy');
for i=1:mm
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:n_x
      fprintf(fid,form1i3g,j,P(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_y - 3*mm) > 0) 
   fprintf(fid,form3i, [3*mm+1:n_y]); fprintf(fid,' := ');
   for j=1:n_x
      fprintf(fid,form1i3g,j,P(j,3*mm+1:n_y));
   end % for
end % if
fprintf(fid,' ;\n');

% ... quadratic part of objective Pyy
fprintf(fid,'\nparam\t \tPyy');
for i=1:mm
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:n_y
      fprintf(fid,form1i3g,j,P(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_y - 3*mm) > 0) 
   fprintf(fid,form3i, [3*mm+1:n_y]); fprintf(fid,' := ');
   for j=1:n_y
      fprintf(fid,form1i3g,j,P(j,3*mm+1:n_y));
   end % for
end % if
fprintf(fid,' ;\n');

% ... linear constraints a, A
fprintf(fid,'\nparam\t: \ta :=');
for i=1:m_1
   fprintf(fid,form1i1g, i, a(i));
end % for
fprintf(fid,' ;\n');
fprintf(fid,'\nparam\t \tAx');
for i=1:nn
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:m_1
      fprintf(fid,form1i3g,j,A(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_x - 3*nn) > 0) 
   fprintf(fid,form3i, [3*nn+1:n_x]); fprintf(fid,' := ');
   for j=1:m_1
      fprintf(fid,form1i3g,j,A(j,3*nn+1:n_x));
   end % for
end % if
fprintf(fid,' ;\n');
if ( norm( A(:,n_x+1:n_x+n_y) ) > 0 )
   fprintf(fid,'\nparam\t \tAy');
   for i=1:mm
      ii = (i-1)*3;
      fprintf(fid,form3ie, [ii+1:ii+3]);
      for j=1:m_1
         fprintf(fid,form1i3g,j,A(j,ii+1:ii+3));
      end % for
      fprintf(fid,'\n');
   end % for
   if ((n_y - 3*mm) > 0) 
      fprintf(fid,form3i, [3*mm+1:n_y]); fprintf(fid,' := ');
      for j=1:m_1
         fprintf(fid,form1i3g,j,A(j,3*mm+1:n_y));
      end % for
   end % if
   fprintf(fid,' ;\n');
end % if

% ... LCP data N*x + M*y + q >= 0
fprintf(fid,'\nparam\t: \tq :=');
for i=1:n_y
   fprintf(fid,form1i1g, i, q(i));
end % for
fprintf(fid,' ;\n');
fprintf(fid,'\nparam\t \tN');
for i=1:nn
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:n_y
      fprintf(fid,form1i3g,j,N(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_x - 3*nn) > 0) 
   fprintf(fid,form3i, [3*nn+1:n_x]); fprintf(fid,' := ');
   for j=1:n_y
      fprintf(fid,form1i3g,j,N(j,3*nn+1:n_x));
   end % for
end % if
fprintf(fid,' ;\n');
fprintf(fid,'\nparam\t \tM');
for i=1:mm
   ii = (i-1)*3;
   fprintf(fid,form3ie, [ii+1:ii+3]);
   for j=1:n_y
      fprintf(fid,form1i3g,j,M(j,ii+1:ii+3));
   end % for
   fprintf(fid,'\n');
end % for
if ((n_y - 3*mm) > 0) 
   fprintf(fid,form3i, [3*mm+1:n_y]); fprintf(fid,' := ');
   for j=1:n_y
      fprintf(fid,form1i3g,j,M(j,3*mm+1:n_y));
   end % for
end % if
fprintf(fid,' ;\n');