function [find, vaf, outperm] = ms_arobfnd(prox, inperm, kblock)

% AROBFND finds and fits an anti-Robinson
% matrix using iterative projection to
% a symmetric proximity matrix in the $L_{2}$-norm based on a 
% permutation identified through the use of iterative quadratic
% assignment.
%
% syntax: [find, vaf, outperm] = arobfnd(prox, inperm, kblock)
%
% PROX is the input proximity matrix ($n \times n$ with a zero main 
% diagonal and a dissimilarity interpretation);
% INPERM is a given starting permutation of the first $n$ integers;
% FIND is the least-squares optimal matrix (with 
% variance-accounted-for of VAF) to PROX having an anti-Robinson 
% form for the row and column object ordering given by the ending
% permutation OUTPERM. KBLOCK defines the block size in the use of the
% iterative quadratic assignment routine.

%tic;
n = size(prox,1);
%targ = targlin(n);

targ = zeros(n,n);

for i = 1:n-1
   for j = (i+1):n
      
      targ(i,j) = abs(i-j);
      targ(j,i) = targ(i,j);
      
   end
end

%[outperm,rawindex,allperms,index] =  order(prox,targlinear,inperm,kblock);

%%%%%%%%%%%%%%%%%%%%%%
n = size(prox,1);
outperm = inperm;
index = 1;

begindex = sum(sum(prox(inperm,inperm).*targ));

for iterate = 1:2
   
   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;
           
         end
         
      end
   end
end

rawindex = begindex;
nchange = 1;

while (nchange == 1)
   
   nchange=0;
   
   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
            end
               
            if ((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
            end
               
           
            
            tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
            if(tryindex > (begindex + 1.0e-008))
               nchange = 1;
               begindex = tryindex;
               outperm = intrperm;
               index = index +1;
               
            end
            
         end
      end
   end
end


rawindex = begindex;
nchange = 1;



while (nchange == 1)
   
   nchange=0;
   
   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;
           
         end
         
      end
   end
end



rawindex = begindex;
nchange = 1;


end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%[fit,vaf] = arobfit(prox,outperm);

inperm = outperm;

%%%%%%%%%%%%%%%%%%%%%
n = size(prox,1);
work = zeros(n*(n-1)*(n-2),1);
fit = prox(inperm,inperm);
cr = 1.0;
proxave = zeros(n,n);

while (cr >= 1.0e-006)
   
   cr = 0.0;
   indexll = 0;
   
   for jone = 1:(n-2)
      for jtwo = (jone+1):(n-1)
         
         p1 = fit(jone,jtwo);
         p2 = fit(jone,jtwo+1);
         
         fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
         fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(indexll+2);
         
         if(fit(jone,jtwo) <= fit(jone,jtwo+1))
            
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
         end
            
         if (fit(jone,jtwo) > fit(jone,jtwo+1))
                        
         ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0;
         work(indexll+1) = ave - fit(jone,jtwo);
         work(indexll+2) = ave - fit(jone,jtwo+1);
         
         fit(jone,jtwo) = ave;
         fit(jone,jtwo+1) = ave;
         
         end
      
         cr = cr + abs(p1-fit(jone,jtwo)) + ...
            abs(p2-fit(jone,jtwo+1));
         
         indexll = indexll + 2;
      end
   end
   
   for jone = 3:n
      for jtwo = 1:(jone-2);
         
         p1 = fit(jtwo,jone);
         p2 = fit(jtwo+1,jone);
         
         fit(jtwo,jone) = fit(jtwo,jone) - work(indexll+1);
         fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(indexll+2);
         
         if(fit(jtwo+1,jone) <= fit(jtwo,jone))
            
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
            
         end
            
         if (fit(jtwo+1,jone) > fit(jtwo,jone))            
         
         ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0;
         work(indexll+1) = ave - fit(jtwo,jone);
         work(indexll+2) = ave - fit(jtwo+1,jone);
         
         fit(jtwo,jone) = ave;
         fit(jtwo+1,jone) = ave;
         
         end
       
         
         cr = cr + abs(p1-fit(jtwo,jone)) + ...
            abs(p2-fit(jtwo+1,jone));
         
         indexll = indexll + 2;
      end
   end
end


for jone = 1:(n-1)
   for jtwo = (jone+1):n
                  
         fit(jtwo,jone) = fit(jone,jtwo);
         
   end
end
   
  aveprox = sum(sum(prox))/(n*(n-1));

for i = 1:n
   for j = 1:n
      if( i ~= j)
         proxave(i,j) = aveprox;
      else
         proxave(i,j) = 0;
      end
   end
end

diff = sum(sum((prox(inperm,inperm) - fit).^2));

denom = sum(sum((prox(inperm,inperm) - proxave).^2));

vaf = 1 - (diff/denom);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nprevperm = 1;

while (nprevperm == 1)
    
    nprevperm = 0;
    prevperm = outperm;
    inperm = outperm;
    
%[outperm,rawindex,allperms,index] = order(prox,fit,inperm,kblock);
  
%%%%%%%%%%%%%%%%
n = size(prox,1);
outperm = inperm;
index = 1;
targ = fit;

begindex = sum(sum(prox(inperm,inperm).*targ));

for iterate = 1:2
   
   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;
           
         end
         
      end
   end
end

rawindex = begindex;
nchange = 1;

while (nchange == 1)
   
   nchange=0;
   
   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
            end
               
            if ((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
            end
               
           
            
            tryindex = sum(sum(prox(intrperm,intrperm).*targ));
            
            if(tryindex > (begindex + 1.0e-008))
               nchange = 1;
               begindex = tryindex;
               outperm = intrperm;
               index = index +1;
               
            end
            
         end
      end
   end
end


rawindex = begindex;
nchange = 1;



while (nchange == 1)
   
   nchange=0;
   
   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;
           
         end
         
      end
   end
end



rawindex = begindex;
nchange = 1;


end    
%%%%%%%%%%%%%%%%%%


%[fit,vaf] = arobfit(prox,outperm);
    
    inperm = outperm;
    
%%%%%%%%%%%%%%%%%%%
n = size(prox,1);
work = zeros(n*(n-1)*(n-2),1);
fit = prox(inperm,inperm);
cr = 1.0;
proxave = zeros(n,n);

while (cr >= 1.0e-006)
   
   cr = 0.0;
   indexll = 0;
   
   for jone = 1:(n-2)
      for jtwo = (jone+1):(n-1)
         
         p1 = fit(jone,jtwo);
         p2 = fit(jone,jtwo+1);
         
         fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
         fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(indexll+2);
         
         if(fit(jone,jtwo) <= fit(jone,jtwo+1))
            
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
         end
            
         if (fit(jone,jtwo) > fit(jone,jtwo+1))
                        
         ave= (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0;
         work(indexll+1) = ave - fit(jone,jtwo);
         work(indexll+2) = ave - fit(jone,jtwo+1);
         
         fit(jone,jtwo) = ave;
         fit(jone,jtwo+1) = ave;
         
         end
      
         cr = cr + abs(p1-fit(jone,jtwo)) + ...
            abs(p2-fit(jone,jtwo+1));
         
         indexll = indexll + 2;
      end
   end
   
   for jone = 3:n
      for jtwo = 1:(jone-2);
         
         p1 = fit(jtwo,jone);
         p2 = fit(jtwo+1,jone);
         
         fit(jtwo,jone) = fit(jtwo,jone) - work(indexll+1);
         fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(indexll+2);
         
         if(fit(jtwo+1,jone) <= fit(jtwo,jone))
            
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
            
         end
            
         if (fit(jtwo+1,jone) > fit(jtwo,jone))            
         
         ave = (fit(jtwo,jone) + fit(jtwo+1,jone))/2.0;
         work(indexll+1) = ave - fit(jtwo,jone);
         work(indexll+2) = ave - fit(jtwo+1,jone);
         
         fit(jtwo,jone) = ave;
         fit(jtwo+1,jone) = ave;
         
         end
       
         
         cr = cr + abs(p1-fit(jtwo,jone)) + ...
            abs(p2-fit(jtwo+1,jone));
         
         indexll = indexll + 2;
      end
   end
end


for jone = 1:(n-1)
   for jtwo = (jone+1):n
                  
         fit(jtwo,jone) = fit(jone,jtwo);
         
   end
end
   
  aveprox = sum(sum(prox))/(n*(n-1));

for i = 1:n
   for j = 1:n
      if( i ~= j)
         proxave(i,j) = aveprox;
      else
         proxave(i,j) = 0;
      end
   end
end

diff = sum(sum((prox(inperm,inperm) - fit).^2));

denom = sum(sum((prox(inperm,inperm) - proxave).^2));

vaf = 1 - (diff/denom);

%%%%%%%%%%%%%%%%%%%
    
    if (sum(prevperm - outperm) ~= 0)
        nprevperm = 1;
    end
end
 
find = fit;

%toc;








