function [fit,vaf,rowperm,colperm] = arobfittm(proxtm,inperm)

%  AROBFITTM does a confirmatory two-mode anti-Robinson fitting of a
%  given ordering of the row and column objects of a two-mode
%  proximity matrix PROXTM using Dykstra's (Kaczmarz's)
%  iterative projection least-squares method.
%
%  syntax: [fit,vaf,rowperm,colperm] = arobfittm(proxtm,inperm)
%
%  INPERM is the given ordering of the row and column objects
%  together; FIT is an nrow (number of rows) by ncol (number of
%  columns) matrix fitted to PROXTM(ROWPERM,COLPERM) 
%  with VAF being the variance-accounted for and
%  based on the (least-squares criterion) sum of 
%  squared discrepancies between FIT and PROXTM(ROWPERM,COLMEAN);
%  ROWPERM and COLPERM are the row and column object orderings
%  derived from INPERM. 

[nrow ncol] = size(proxtm);
n = nrow + ncol;
work = zeros(n*(n-1)*(n-2)*(n-3),1);
tempfit = zeros(n,n);
fit = zeros(nrow,ncol);

for irow = 1:nrow
   for icol = 1:ncol
      
      tempfit(irow,icol+nrow) = proxtm(irow,icol);
      tempfit(icol+nrow,irow) = proxtm(irow,icol);
      
   end
end

mfit = tempfit(inperm,inperm);

cr = 1.0;

while(cr >= 1.0e-005)
   
   cr = 0.0;
   indexll=0;
   
   for ione = 1:(n-2)
      for itwo = (ione+1):(n-1)
         for ithree = (itwo+1):n
                          
               if((((inperm(ione) <= nrow) && (inperm(itwo) <= nrow)) ...
                     && ((inperm(ithree) > nrow) )) || ...
                     (((inperm(ione) > nrow) && (inperm(itwo) > nrow)) ...
                     && ((inperm(ithree) <= nrow) )))
                  
                  p1 = mfit(ione,ithree);
                  p2 = mfit(itwo,ithree);
                                    
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+1);
                  mfit(itwo,ithree) = mfit(itwo,ithree) - work(indexll+2);
                 
                  if (mfit(ione,ithree) >= mfit(itwo,ithree))
                      del = 0;
                  else
                      del = (mfit(itwo,ithree) - mfit(ione,ithree))/2;
                  end
                  
                               
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(itwo,ithree) = mfit(itwo,ithree) - del;
                  
                  
                  work(indexll+1) = del;
                  work(indexll+2) = -del;
                  
                  
                  cr = cr + abs(p1 - mfit(ione,ithree)) + ...
                     abs(p2 - mfit(itwo,ithree));
                     
                 indexll = indexll + 2;
                  
                
                elseif((((inperm(ione) <= nrow) && (inperm(itwo) > nrow)) ...
                  && ((inperm(ithree) > nrow) )) || ...
                  (((inperm(ione) > nrow) && (inperm(itwo) <= nrow)) ...
                  && ((inperm(ithree) <= nrow) )))
                  
                  p1 = mfit(ione,itwo);
                  p2 = mfit(ione,ithree);
                  
                  
                  mfit(ione,itwo) = mfit(ione,itwo) - work(indexll+1);
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+2);
                  
                  if (mfit(ione,ithree) >= mfit(ione,itwo))
                      del = 0;
                  else
                      del = (mfit(ione,itwo) - mfit(ione,ithree))/2;
                  end
                  
                               
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(ione,itwo) = mfit(ione,itwo) - del;
                  
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = del;
                  
                  
                  cr = cr + abs(p1 - mfit(ione,itwo)) + ...
                     abs(p2 - mfit(ione,ithree));
                     
                 indexll = indexll + 2;
                  
                  
                  
               end
               
            end
         end
      end
    
   
end

   
   
   for i = 1:(n-1)
      for j = (i+1):n
         
         mfit(j,i) = mfit(i,j);
         
      end
   end
   
   rowidx = 0;
   
   for i = 1:n
      
      if(inperm(i) <= nrow)
         
         rowidx = rowidx + 1;
         
         colidx = 0;
         
         for j = 1:n
            
            if(inperm(j) > nrow)
               
               colidx = colidx + 1;
               fit(rowidx,colidx) = mfit(i,j);
            end
         end
      end
   end
   
   rowperm = zeros(nrow,1);
   colperm = zeros(ncol,1);
   rowidx = 0;
   colidx = 0;
   
   for i = 1:n
      
      if(inperm(i) <= nrow)
         
         rowidx = rowidx + 1;
         rowperm(rowidx) = inperm(i);
         
      elseif(inperm(i) > nrow)
         
         colidx = colidx + 1;
         colperm(colidx) = inperm(i) - nrow;
         
      end
   end
   
   aveprox = sum(sum(proxtm))/(nrow*ncol);
   
   diff = sum(sum((proxtm(rowperm,colperm) - fit).^2));
   denom = sum(sum((proxtm(rowperm,colperm) - aveprox).^2));
   
   vaf = 1 - (diff/denom);
  
  
   
       