function [find,vaf,targone,targtwo,outpermone,outpermtwo, ...
        rowpermone,colpermone,rowpermtwo,colpermtwo] = ...
    biarobfndtm(proxtm,inpermone,inpermtwo,kblock)

% BIAROBFNDTM finds and fits the sum of two anti-Robinson using iterative projection to
% a two-mode proximity matrix in the $L_{2}$-norm based on permutations 
% identified through the use of iterative quadratic assignment.
% PROXTM is the input two-mode proximity matrix ($nrow \times ncol$)
% and a dissimilarity interpretation);
% FIND is the least-squares optimal matrix (with variance-accounted-for
% of VAF) to PROXTM and is the sum of the two  matrices
% TARGONE and TARGTWO based on the two row and column
% object orderings given by the ending permutations OUTPERMONE 
% and OUTPERMTWO, and in turn ROWPERMONE and ROWPERMTWO and COLPERMONE
% and COLPERMTWO. KBLOCK defines the block size in the use the
% iterative quadratic assignment routine; the input permutations are INPERMONE and
% INPERMTWO.  

[nrow ncol] = size(proxtm);
n = nrow + ncol;


[targ1,vaftarg1,outperm1,rowperm1,colperm1] = arobfndtm(proxtm,inpermone,kblock);

for i = 1:nrow
    for j = 1:ncol
        
          resprox1(rowperm1(i),colperm1(j)) = proxtm(rowperm1(i),colperm1(j)) ...
             - targ1(i,j);
     
     
     
 end
end
            
[targ2,vaftarg2,outperm2,rowperm2,colperm2] = arobfndtm(resprox1,inpermtwo,kblock);

for i = 1:nrow
    for j = 1:ncol
        
          resprox2(rowperm2(i),colperm2(j)) = resprox1(rowperm2(i),colperm2(j)) ...
             - targ2(i,j);
     
     
     
 end
end

find = proxtm - resprox2;


aveprox = sum(sum(proxtm))/(nrow*ncol);

for i = 1:nrow
    for j = 1:ncol
        
        proxave(i,j) = aveprox;
        
    end
    
end

diff = sum(sum((proxtm - find).^2));
denom = sum(sum((proxtm - proxave).^2));
vaf = 1 - (diff/denom);
            
targone = targ1;
targtwo = targ2;
rowpermone = rowperm1;
rowpermtwo = rowperm2;
colpermone = colperm1;
colpermtwo = colperm2;
outpermone = outperm1;
outpermtwo = outperm2;

vafdiff = 1.0;
iteration = 0;

while ((vafdiff >= 1.0e-004) & (iteration < 20))
    
    vafprev = vaf;
    iteration = iteration + 1;
    
   for i = 1:nrow
    for j = 1:ncol
        
           resprox(rowpermtwo(i),colpermtwo(j)) = proxtm(rowpermtwo(i),colpermtwo(j)) ...
             - targtwo(i,j);
          
 end
end

          [targone,vafone,outpermone,rowpermone,colpermone] = arobfndtm(resprox,inpermone,kblock);
    
  for i = 1:nrow
    for j = 1:ncol
       
        resprox(rowpermone(i),colpermone(j)) = proxtm(rowpermone(i),colpermone(j)) ...
             - targone(i,j);
     
    
     
 end
end

      [targtwo,vaftwo,outpermtwo,rowpermtwo,colpermtwo] = arobfndtm(resprox,inpermtwo,kblock);

  for i = 1:nrow
    for j = 1:ncol
        
        find(rowpermone(i),colpermone(j)) = targone(i,j);
     
    
 end
end

for i = 1:nrow
    for j = 1:ncol
        
            find(rowpermtwo(i),colpermtwo(j)) = find(rowpermtwo(i),colpermtwo(j)) ...
             + targtwo(i,j);
     
     
     
 end
end

   
    diff = sum(sum((proxtm - find).^2));
    denom = sum(sum((proxtm - proxave).^2));
    vaf = 1 - (diff/denom);
    
        
    vafdiff = abs(vaf - vafprev);
    
end


