function [fit,diff,coord,exitflag] = ...
   linfitl1(prox,inperm)

%LINFITL1 does a confimatory fitting in the $L_{1}$ norm of a given unidimensional
%  order using linear programming.
%  INPERM is the given order; FIT is an $n \times n$ matrix that is fitted to
%  PROX(INPERM,INPERM) with $L_{1}$ value DIFF; COORD gives the ordered coordinates
%  whose absolute differences could be used to reconstruct FIT.
%  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).


n=size(prox,1);
coord = zeros(n,1);
proxmat = prox(inperm,inperm);
nch2 = n*(n-1)/2;
xlength = 3*nch2 + n;
aequal = zeros(2*nch2+1,xlength);

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

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

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

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

bconst(2*nch2+1) = 0;

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

fweight = zeros(xlength,1);

for i = 1:xlength
   if((i > n) & (i < (2*nch2+n+1)))
      fweight(i) = 1;
   end
end

options = optimset('MaxIter',100000);
[xend,diff,exitflag,output] = ...
linprog(fweight,[],[],aequal,bconst,lbound,[],[],options);

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

dev = 0;
fit = zeros(n,n);
for i = 1:n
    for j = 1:n
        if (i ~= j)
            fit(i,j) = abs(coord(i) - coord(j));
            dev = dev + abs(proxmat(i,j) - fit(i,j));
        end
    end
end

if (exitflag <= 0)
    exitflag
end

