function [startcoord, begval, outcoord, endval, exitflag] = ...
   unifitl2nlp(prox,inperm)

%UNIFITL2NLP carries out a unidimensional scaling of a symmetric proximity 
%  matrix using the l2 norm and the nonlinear programming refomulation of 
%  Lau, Leung, and Tse.
%  PROX is the input proximity matrix (with a zero main diagonal and a
%  dissimilarity interpretation);
%  INPERM is the input permutation of the first $n$ integers which is used to
%  obtain the starting coordinates (STARTCOORD) and the initial least squares 
%  loss value (BEGVAL);
%  OUTCOORD are the ending coordinates values having the final least squares loss
%  value of ENDVAL;
%  EXITFLAG indicates the success of the optimization ( > 0 indicates convergence;
%  0 indicates that the maximum number of function evaluations or iterations were 
%  reached; and < 0 denotes nonconvergence).
%
%  syntax: [startcoord, begval, outcoord, endval, exitflag] = ...
%   unifitl2nlp(prox,inperm)

n = size(prox,1);
outcoord = zeros(n,1);
tmpcoord = zeros(n,1);
coord = zeros(n,1);
proxmat = prox(inperm,inperm);

for i = 1:n
   for j = 1:n
      if(i>j)
         tmpcoord(i) = tmpcoord(i) + proxmat(i,j);
      end
      if(i<j)
         tmpcoord(i) = tmpcoord(i) - proxmat(i,j);
      end
   end
   
   tmpcoord(i) = tmpcoord(i)/n;
   
end

for i = 1:n
    for j = 1:n
        if(inperm(i) == j)
            coord(j) = tmpcoord(i);
        end
    end
end

nch2 = n*(n-1)/2;
xlength = (4*nch2) + n;
xstart = zeros(xlength,1);


for i = 1:n
   xstart(i) = coord(i);
end

index = n;
for i = 2:n
   for j = 1:(i-1)
      index = index + 1;
      xstart(index) = prox(i,j) - coord(i) + coord(j);
      xstart(index+nch2) = prox(i,j) - coord(j) + coord(i);
      if(coord(i) > coord(j))
         xstart(index+2*nch2) = 1;
      end
      if(coord(i) < coord(j))
         xstart(index+3*nch2) = 1;
      end
   end
end


aequal = zeros(3*nch2+1,xlength);

index = 0;
for i = 2:n
   for j = 1:(i-1)
      index = index + 1;
      aequal(index,i) = 1;
      aequal(index,j) = -1;
      aequal(index,index+n) = 1;
      aequal(index+nch2,i) = -1;
      aequal(index+nch2,j) = 1;
      aequal(index+nch2,index+nch2+n) = 1;
      aequal(index+2*nch2,index+2*nch2+n) = 1;
      aequal(index+2*nch2,index+3*nch2+n) = 1;
   end
end

for i = 1:n
   aequal(3*nch2+1,i) = 1;
end


bconst = zeros(3*nch2+1,1);

index = 0;
for i = 2:n
   for j = 1:(i-1)
      index = index + 1;
      bconst(index) = prox(i,j);
      bconst(index+nch2) = prox(i,j);
      bconst(index+2*nch2) = 1;
   end
end

bconst(3*nch2+1) = 0;

for i = 1:xlength
   lbound(i) = -inf;
   if(i > (n+2*nch2))
      lbound(i) = 0;
   end
end


[xend,fval,exitflag,output] = ...
   fmincon('objfunl2',xstart,[],[],aequal,bconst,lbound,[]);
for i = 1:n
   outcoord(i) = xend(i);
end


startcoord = coord;
begval = objfunl2(xstart);
outcoord;
endval = objfunl2(xend);

