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


%UNIFITL1NLP carries out a unidimensional scaling of a symmetric proximity 
%  matrix using the $L_{1}$ 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 $L_{1}$ norm 
%  loss value (BEGVAL);
%  OUTCOORD are the ending coordinates values having the final $L_{1}$ 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).


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

   
[fit,diff,tmpcoord,exitflag] = linfitl1(proxmat,1:n)

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

nch2 = n*(n-1)/2;
xlength = (6*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;
      e1 = prox(i,j) -coord(i) + coord(j);
      if(e1 >= 0)
         xstart(index) = e1;
         xstart(index+nch2) = 0;
      end
       
      if(e1 < 0)
         xstart(index) = 0;
         xstart(index+nch2) = -e1;
      end
      e2 = prox(i,j) -coord(j) + coord(i);
      if(e2 >= 0)
         xstart(index+2*nch2) = e2;
         xstart(index+3*nch2) = 0;
      end
      if(e2 < 0)
         xstart(index+2*nch2) = 0;
         xstart(index+3*nch2) = -e2;
      end
      
      if(coord(i) >= coord(j))
         xstart(index+4*nch2) = 1;
      end
      if(coord(i) < coord(j))
         xstart(index+5*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,index+n+nch2) = -1;
      aequal(index+nch2,i) = -1;
      aequal(index+nch2,j) = 1;
      aequal(index+nch2,index+2*nch2+n) = 1;
      aequal(index+nch2,index+3*nch2+n) = -1;
      aequal(index+2*nch2,index+4*nch2+n) = 1;
      aequal(index+2*nch2,index+5*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)
      lbound(i) = 0;
   end
end



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

startcoord = coord;
begval = objfunl1(xstart);
endval = objfunl1(xend);

toc