%
% [fn, numlabelconst, totalnumconstraints] = printSDPAord(
%    fid, y, maxoravg, C, perrowthresh, requirethreshord, comment)
%
% Sets up the SDP for max-margin matrix factorization with ordinal
% labels and prints the problem in sparse SDPA format to a file.  This
% is the same problem as solveDord, but instead of solving it with
% YALMIP, it is just printed to a file.  You should then use a SDP
% solver on this problem, and read the solution using readSDPAord.
%
% Inputs:
%
% fid = either a file ID to which the problem is printed, or a
% string, which is used as a base for a filename.  If it is a
% string, a filename is created by appending to it suffixes
% describing the parameters, and the full filename is returned in
% the output argument.
%
% y = a matrix of labels.  Labels are (small) positive integers
% 1..R (the range is inferred from y).  0 indicates a missing label.
%
% maxoravg = 'a' for nuclear-norm, 'm' for max-norm (default is
% nuclear-norm)
% 
% C = label loss function
%   C>0: Shashua-Levin type hinge loss on immediate threshold, scaled by C
%   C<0: Hinge loss summed over all thresholds, scaled by (-C)
%   C=inf and C=-inf force hard margin constraints, but lead to a
%   different SDP
%   Default is C=inf, i.e. no slack, with constraints only on immediate
%   thresholds
%
% perrowthresh = 1 to allow different thresholds for each row. Default is
% 0, i.e. universal thresholds for entire matrix
% 
% requirethreshord = 1 to require thresholds be correctly ordered.
% Changing this is highly unrecommended and might lead to very strange
% SDPs.
%
% comment = an optional string written as a comment in the header.
%
% Copyright May 2004, Nathan Srebro, nati@mit.edu
%

function [fn, numlabelconst, totalnumconstraints] = printSDPAord(...
    fid, y, maxoravg, C, perrowthresh, requirethreshord, comment)
  
  [n,m] = size(y);
  [i,a,v] = find(y);
  p=length(v);
  R=max(v);
  if (nargin>2) & (maxoravg(1)=='m')
    maxoravg = 'max';
    maxprob = 1;
  else
    maxoravg = 'avg';
    maxprob = 0;
  end
  if (nargin<4)
    C = inf;
  end
  if C<0
    C = -C;
    sumrankmarg = 1;
  else
    sumrankmarg = 0;
  end
  if nargin<5
    perrowthresh = 0;
  end
  if nargin<6
    requirethreshord = 1;
    % It is highly unrecomneded to set this otherwise, unless allth is
    % used and there is a label between every two thresholds
  end
  if sumrankmarg == 1
    fnpenaltytype = 'allth';
    compenaltytype = 'all thresholds';
  else
    fnpenaltytype = 'imdth';	% imidiate thresholds, ordered 
    compenaltytype = 'imidiate thresholds';
  end
  allowslack = (C<inf);
  if perrowthresh
    perrowthreshtext = 'rowth';
    numthvecs = n;
  else
    perrowthreshtext = 'unith';
    numthvecs = 1;
  end
  if requirethreshord
    fnreqthordtext = '';
    comreqthordtext = '';
  else
    fnreqthordtext = 'UNC';
    comreqthordtext = 'ord UNCONSTRAINED';
  end

  if ischar(fid)
    fn = sprintf('%s.%s_%s_%s_%s%s.dat-s',fid,maxoravg,fnpenaltytype,num2str(C),perrowthreshtext,fnreqthordtext);
    fid = fopen(fn,'w');
    remembertoclose = 1;
  else
    fn = [];
    remembertoclose = 0;
  end
  
  % Comments
  fprintf(fid,'* MMMF-MC, %dx%d, %d binnary labels (%d sparse)\n',n,m,p, ...
	  p/n/m);
  fprintf(fid,'* Minimizing %s-norm, %s hinge loss penalty %f, %s %s\n',...
	  maxoravg, compenaltytype, C, perrowthreshtext, comreqthordtext);
  fprintf(fid,'* Generated by printSDPA3, %s\n',date);
  fprintf(fid,'* printSDPA3: Nathan Srebro, May 2004, nati@mit.edu\n');
  if nargin>6
    fprintf(fid,'* %s\n',comment);
  end

  % bound on thresholds: all thresholds will be in -thbound<...<thbound
  thbound = R*10;
  
  % How many constraints?
  % Label constraints
  if sumrankmarg
    numlabelconst = (R-1)*p;
  else
    numlabelconst = 2*p-(sum(v==1)+sum(v==R));
  end
  % Max-norm constraints
  normconstofset = numlabelconst;
  nummaxprobconst = maxprob*(n+m-1);
  % Threshold order constraints (also bound on highest threshold)
  thconstofset = nummaxprobconst+normconstofset;
  numthordconst = requirethreshord*(R-2)*numthvecs;
  numthboundconst = numthvecs;
  numthconst = numthordconst+numthboundconst;
  % Bound overall negative bias
  biasboundconstofset = thconstofset+numthconst;
  numbiasboundconst = 1;
  
  % Number of primal constraints / dual variables
  totalnumconstraints = numlabelconst+nummaxprobconst+numthconst+numbiasboundconst;  
  fprintf(fid,'%d\n',totalnumconstraints); 

  % The blocks of each primal constrain matrix / dual var matrix:
  % (can be interpeted as ndual matrix constraints, or number of
  % primal matrix variables or variable groups)

  % First block correspond to [ () X ; X' () ]
  % Second block is a single overall negative bias
  % Third block is thresholds
  % Fourth block is positive distance of points from margin
  blocksizes = [(n+m) , -1, -numthvecs*(R-1) , -numlabelconst];
  numblocks=4;
  % If slack is alowed, Fifth block is slack, ie negative distance from margin
  if allowslack
    numblocks=numblocks+1;
    slackblock = numblocks;
    blocksizes = [blocksizes , -numlabelconst];
  end
  % If threshold order is constrained, next block is distance between
  % threholds.  Otherwise, it is just used to bound the highest threshold
  numblocks=numblocks+1;
  thordblock = numblocks;
  blocksizes = [blocksizes , -numthconst];

  
  % Number of blocks
  fprintf(fid,'%d\n',numblocks);
  % The size of each such block: (negative indicated diagonal block)
  fprintf(fid,'%d ',blocksizes);
  fprintf(fid,'\n');

  % Primal free terms in constraints / dual objective coeficients
  fprintf(fid,'%f ',[ones(1,numlabelconst) zeros(1,nummaxprobconst+numthordconst) ...
		     repmat(2*thbound,1,numthboundconst) thbound]);
  fprintf(fid,'\n');

  % Now come the matrices, formated <mat#> <blk#> <i> <j> <value>
  % These are the primal constraints, and the dual mat struct

  % Set overall negative bias term
  fprintf(fid,'%d 2 1 1 1\n',biasboundconstofset+1);
  
  % label constraints -- terms dependent on labels (X_ia and threholds)
  rowthstride = perrowthresh*(R-1);
  if sumrankmarg
    vv = sign(2*op(v,'>',1:(R-1))-1) ;
    % X_ia terms    
    fprintf(fid,'%d 1 %d %d %f\n',[1:(p*(R-1)) ; repmat(i',1,R-1) ; repmat(n+a',1,R-1) ; ...
		    vv(:)'/2]);
    % nagative bias
    fprintf(fid,'%d 2 1 1 %f\n',[1:(p*(R-1)) ; vv(:)']);
    % thresholds
    whichth = op(1:(R-1),'+',rowthstride*(i-1)) ;
    fprintf(fid,'%d 3 %d %d %f\n',[1:(p*(R-1)) ; whichth(:)' ; ...
		    whichth(:)' ; -vv(:)']) ;
  else
    curconst = 0;
    for rr=1:(R-1)
      left = (v==rr);
      right = (v==(rr+1));
      % X_ia term for labels imidiately to the left
      fprintf(fid,'%d 1 %d %d -0.5\n',[curconst+(1:nnz(left)) ; i(left)' ; ...
		    n+a(left)' ]);
      % nagative bias
      fprintf(fid,'%d 2 1 1 -1\n',curconst+(1:nnz(left)));
      % thresholds for them
      fprintf(fid,'%d 3 %d %d 1.0\n',[curconst+(1:nnz(left)) ; rr+ ...
		    rowthstride*(i(left)-1)' ; rr+rowthstride*(i(left)-1)' ]) ;
      curconst = curconst + nnz(left);
      % X_ia term for labels imidiately to the right
      fprintf(fid,'%d 1 %d %d 0.5\n',[curconst+(1:nnz(right)) ; i(right)' ; ...
		    n+a(right)' ]);
      % nagative bias
      fprintf(fid,'%d 2 1 1 1\n',curconst+(1:nnz(right)));
      % thresholds for them
      fprintf(fid,'%d 3 %d %d -1.0\n',[curconst+(1:nnz(right)) ; rr+ ...
		    rowthstride*(i(right)-1)' ; rr+rowthstride*(i(right)-1)' ]) ;
      curconst = curconst + nnz(right);
    end
    if numlabelconst ~= curconst
      disp('ERROR!!! Missmatch in constraint counting!');
      numlabelconst, curconst
    end
  end

  % positive distance from margin terms
  fprintf(fid,'%d 4 %d %d -1\n',repmat(1:numlabelconst,3,1));
  % slack (negative distance from margin)
  if allowslack
    fprintf(fid,'%d 5 %d %d 1\n',repmat(1:numlabelconst,3,1));
    % Add slack to objective / bound dual vars by C
    fprintf(fid,'0 5 %d %d %f\n',[1:numlabelconst ; 1:numlabelconst ; repmat(-C,1,numlabelconst)]);
  end
  
  if maxprob   % Max-norm constraint
    % In primal: all norms are equal (bounded and equal are equiv for opt)
    % this is represented by requiring all elements on diagonal of first
    % block to be equal to the leading element on the diagonal.
    % In the dual, this adds variables along the diagonal
    fprintf(fid,'%d 1 %d %d 1.0\n', ...
	    [normconstofset+(1:(n+m-1)); repmat(2:(n+m),2,1)]);
    fprintf(fid,'%d 1 1 1 -1.0\n', normconstofset+(1:(n+m-1)));
    % Primal: Objective is (negative) leading element on diagonal
    fprintf(fid,'0 1 1 1 -1.0\n');
    % In dual: trace of 1st block is equal to one
  else % Avarage-norm constraint
    % objective is sum of norms (in primal)
    % diagonal of first block is all ones (in dual)
    fprintf(fid,'0 1 %d %d -1.0\n',[1:(n+m) ; 1:(n+m)]);
  end

  if requirethreshord			% Threshold order constraints
    negthindx = op(2:(R-1),'+',rowthstride*(0:(numthvecs-1))');
    posthindx = op(1:(R-1),'+',rowthstride*(0:(numthvecs-1))');
  else
    negthindx = [];
    posthindx =  (R-1)+rowthstride*(0:(numthvecs-1))';
  end    
  fprintf(fid,'%d 3 %d %d -1.0\n',[thconstofset+(1:numthordconst)  ; ...
		    negthindx(:)' ; negthindx(:)' ]);
  fprintf(fid,'%d 3 %d %d 1.0\n',[thconstofset+(1:numthconst)  ; ...
		    posthindx(:)' ; posthindx(:)' ]);
  fprintf(fid,'%d %d %d %d 1.0\n',[thconstofset+(1:numthconst)  ; ...	
		    repmat(thordblock,1,numthconst) ; ...
		    1:numthconst ; 1:numthconst ]);
  
  if remembertoclose
    fclose(fid);
  end
  
  
  
function out = op(arg1,operator,arg2)
  % computes the binnary operation on the arguments, extending 1-dim
  % dimmensions apropriately. E.g. it is ok to multiply 1xNxP and
  % MxNx1 matrices, subtarct a vector from a matrix, etc.
  %
  % Written by Nathan Srebro, MIT LCS, October 1998.
  
  shape1=[size(arg1),ones(1,length(size(arg2))-length(size(arg1)))] ;
  shape2=[size(arg2),ones(1,length(size(arg1))-length(size(arg2)))] ;
  
  out = feval(operator, ...
      repmat(arg1,(shape1==1) .* shape2 + (shape1 ~= 1) ), ...
      repmat(arg2,(shape2==1) .* shape1 + (shape2 ~= 1) )) ;