function [outperm, rawindex, allperms, index, coord, diff] = ...
   uniscalqa(prox, targ, inperm, kblock)

%UNISCALQA carries out a unidimensional scaling of a symmetric proximity
%  matrix using iterative quadratic assignment.
%  PROX is the input proximity matrix (with a zero main diagonal and a 
%  dissimilarity interpretation);
%  TARG is the input target matrix (usually with a zero main diagonal and 
%  with a dissimilarity interpretation representing equally-spaced locations 
%  along a continuum);
%  INPERM is the input beginning permutation (a permuation of the first $n$ integers).
%  OUTPERM is the final permutation of PROX with the cross-product index RAWINDEX
%  with respect to TARG redefined as $ = \{abs(coord(i) - coord(j))\}$;
%  ALLPERMS is a cell array containing INDEX entries corresponding to all the
%  permutations identified in the optimization from ALLPERMS{1} = INPERM to 
%  ALLPERMS{INDEX} = OUTPERM.
%  The insertion and rotation routines use from 1 to KBLOCK
%  (which is less than or equal to $n-1$) consecutive objects in 
%  the permutation defining the row and column order of the data matrix.
%  COORD is the set of coordinates of the unidimensional scaling 
%  in ascending order;
%  DIFF is the value of the least-squares loss function for the
%  coordinates and object permutation.

tic;
n = size(prox,1);
outperm = inperm;
index = 1;
allperms{index} = inperm;
begindex = sum(sum(prox(inperm,inperm).*targ));
coord = zeros(n,1);

prevperm = inperm;
nprevperm = 1;

while (nprevperm == 1)
    nprevperm = 0;

    nchange = 1;
   
while (nchange == 1)
   
   nchange=0;
   
   for k = 1:(n-1)
      for j = (k+1):n
                   
         intrperm = outperm;
             
         intrperm(k) = outperm(j);
         intrperm(j) = outperm(k);
         
         tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
         if(tryindex > (begindex + 1.0e-008))
            nchange = 1;
            begindex = tryindex;
            outperm = intrperm;
            index = index + 1;
            allperms{index} = intrperm;
         end
         
      end
   end
   
   
   for k = 1:kblock
      for insertpt = 1:(n+1)
         for nlimlow = 1:(n+1-k)
            
            intrperm = outperm;
            
            if (nlimlow > insertpt)
               
               jtwo = 0;
               for j = insertpt:(insertpt+k-1)
                  intrperm(j) =outperm(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = (insertpt+k):(nlimlow+k-1);
                  intrperm(j) = outperm(insertpt+jone);
                  jone = jone + 1;
               end
               
            elseif ((nlimlow+k) < insertpt)
               
               jtwo = 0;
               for j = (insertpt-k):(insertpt-1)
                  intrperm(j) = outperm(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = nlimlow:(insertpt-k-1)
                  intrperm(j) = outperm(nlimlow+k+jone);
                  jone = jone + 1;
               end
               
            else
               
            end
            
            tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
            if(tryindex > (begindex + 1.0e-008))
               nchange = 1;
               begindex = tryindex;
               outperm = intrperm;
               index = index +1;
               allperms{index} = intrperm;
            end
            
         end
      end
   end
   

if (kblock > 1)
    
   for k = 2:kblock
      for nlimlow = 1:(n+1-k)
                   
         intrperm = outperm;
            
         for j = 1:k
            intrperm(nlimlow+j-1) = outperm(nlimlow+k-j);
         end
         
         tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
         if(tryindex > (begindex + 1.0e-008))
            nchange = 1;
            begindex = tryindex;
            outperm = intrperm;
            index = index + 1;
            allperms{index} = intrperm;
         end
         
      end
   end
end

end





rawindex = begindex;
if(any(prevperm - outperm) == 1)
    nprevperm = 1;
    prevperm = outperm;
end

coord = zeros(n,1);

for i = 1:n
   for j = 1:n
      if(i > j)
         coord(i) = coord(i) + prox(outperm(i),outperm(j));
      end
      if(i < j) 
         coord(i) = coord(i) - prox(outperm(i),outperm(j));
      end
   end
   coord(i) = coord(i)/n;
end

diff = 0;
for i = 1:(n-1)
   for j = (i+1):n
      diff = diff + (prox(outperm(i),outperm(j)) - ...
         abs(coord(i) - coord(j)))^2;
   end
end

for i = 1:n
    for j = 1:n
        targ(i,j) = abs(coord(i) - coord(j));
    end
end
  
begindex = sum(sum(prox(outperm,outperm).*targ));

end

toc
