function [pcoordsort, pperm, gcoordsort, gperm, gdiff, pdiff] = ...
   plinorder(prox,inperm)

%PLINORDER carries out a unidimensional scaling of a symmetric proximity
%  matrix using Pliner smoothing function in Guttman's updating algorithm.
%  PROX is the input proximity matrix (with a zero main diagonal and a
%  dissimilarity interpretation);
%  INPERM is an input permutation of the first $n$ integers;
%  GCOORDSORT are the coordinates from the Guttman update ordered from most negative to most positive;
%  GPERM is the object permutation that indicates which objects are at which of the
%  ordered coordinates in GCOORDSORT;
%  PCOORDSORT are the coordinates from the Pliner smoother of the Guttman update ordered from most
%  negative to most positive;
%  PPERM is the object permutation that indicates which objects are at which of the ordered
%  coordinates in PCOORDSORT;
%  GDIFF is the value of the least-squares loss function for the 
%  coordinates and object permutation obtained from the Guttman update;
%  PDIFF is the value of the least-squares loss function for the coordinates 
%  and object permutation obtained from the Pliner smoothing procedure.
%
% syntax: [pcoordsort, pperm, gcoordsort, gperm, gdiff, pdiff] = plinorder(prox,inperm)

n = size(prox,1);
gpermprev = inperm;
gcoord = zeros(n,1);
gcoordprev = zeros(n,1);

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


for i = 1:n
   for j = 1:n
      gcoord(i) = gcoord(i) + (prox(i,j)*sign(gcoordprev(i) - ...
         gcoordprev(j)))/n;
   end
end

gperm = 1:n;
gcoordsort = gcoord;

for i = 1:(n-1)
   for j = (i+1):n
      if(gcoordsort(i) > gcoordsort(j))
         tempcoord = gcoordsort(i);
         gcoordsort(i) = gcoordsort(j);
         gcoordsort(j) = tempcoord;
         tempperm = gperm(i);
         gperm(i) = gperm(j);
         gperm(j) = tempperm;
      end
   end
end
gpermprev = randperm(n);

while(any(gpermprev - gperm) == 1)
   gcoordprev = gcoord;
   gpermprev = gperm;
   gcoord = zeros(n,1);   
   for i = 1:n
      for j = 1:n
         gcoord(i) = gcoord(i) + (prox(i,j)*sign(gcoordprev(i) - ...
            gcoordprev(j)))/n;
      end
   end
   
  
gperm = 1:n;
gcoordsort = gcoord;

   for i = 1:(n-1)
      for j = (i+1):n
         if(gcoordsort(i) > gcoordsort(j))
            tempcoord = gcoordsort(i);
            gcoordsort(i) = gcoordsort(j);
            gcoordsort(j) = tempcoord;
            tempperm = gperm(i);
            gperm(i) = gperm(j);
            gperm(j) = tempperm;
         end
      end
   end
   
   gdiff = 0;
   for i = 1:n
      for j = 1:n
         gdiff = gdiff + (prox(i,j) - abs(gcoord(i) - gcoord(j)))^2;
      end
   end
  end



ep = 2*(max(sum(prox))/n);
pperm = gperm;
pcoord = gcoord;

for k = 2:100
      
   ppermprev = randperm(n);
   
   while(any(ppermprev - pperm) == 1)

   pcoordprev = pcoord;
   ppermprev = pperm;
   
   pcoord = zeros(n,1);
   for i = 1:n
      for j = 1:n
         abst = abs(pcoordprev(i) - pcoordprev(j));
         if(abst < ep)
            factor = ((pcoordprev(i) - pcoordprev(j))/ep)* ...
               (2 - (abst/ep));
         end
         if(abst >= ep)
             factor = sign(pcoordprev(i) - pcoordprev(j));
         end
            pcoord(i) = pcoord(i) + (prox(i,j)*factor)/n;
      end
   end
   
   
  
pperm = 1:n;
pcoordsort = pcoord;

   for i = 1:(n-1)
      for j = (i+1):n
         if(pcoordsort(i) > pcoordsort(j))
            tempcoord = pcoordsort(i);
            pcoordsort(i) = pcoordsort(j);
            pcoordsort(j) = tempcoord;
            tempperm = pperm(i);
            pperm(i) = pperm(j);
            pperm(j) = tempperm;
         end
      end
   end
   
   
   pdiff = 0;
   for i = 1:n
      for j = 1:n
         pdiff = pdiff + (prox(i,j) - abs(pcoord(i) - pcoord(j)))^2;
      end
   end


end
ep = ep*(100 - k +1)/100;


end
gdiff = .5*gdiff;
pdiff = .5*pdiff;



