function [fit, vaf] = arobfit(prox, inperm)

% AROBFIT fits an anti-Robinson matrix using iterative projection to        
% a symmetric proximity matrix in the $L_{2}$-norm.
%
% syntax: [fit, vaf] = arobfit(prox, inperm)
%
% PROX is the input proximity matrix ($n \times n$ with a zero main
% diagonal and a dissimilarity interpretation);
% INPERM is a given permutation of the first $n$ integers;
% FIT 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 INPERM.

%tic;
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;
            
         elseif(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;
            
         elseif(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);

%toc;         