function [find,vaf,targone,targtwo,outpermone,outpermtwo, ...
        rowpermone,colpermone,rowpermtwo,colpermtwo,addconone,...
        addcontwo,coordone,coordtwo,axes] = ...
    biscaltmac(proxtm,inpermone,inpermtwo,kblock,nopt)

% BISCALTMAC finds and fits the sum of two linear unidimensional scales 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 and ADDCONONE and ADDCONTWO are
% the two additive constants for the two model components; The $n$ coordinates
% are in COORDONE and COORDTWO.  The input permutations are INPERMONE and
% INPERMTWO.  The $n \times 2$ matrix AXES gives the plotting coordinates for the
% combined row and column object set.
%  NOPT controls the confirmatory or exploratory fitting of the unidimensional
%  scales; a value of NOPT = 0 will fit in a confirmatory manner the two scales
%  indicated by INPERMONE and INPERMTWO; a value of NOPT = 1 uses iterative QA
%  to locate the better permutations to fit;

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

if (nopt == 1)
    
      [targ1,vaftarg1,outperm1,rowperm1,colperm1,addcon1,coord1] = uniscaltmac(proxtm,inpermone,kblock);

end

if (nopt == 0)
    
    outperm1 = inpermone;
    [targ1,vaftarg1,rowperm1,colperm1,addcon1,coord1] = linfittmac(proxtm,inpermone);
    
end


for i = 1:nrow
    for j = 1:ncol
        
          resprox1(rowperm1(i),colperm1(j)) = proxtm(rowperm1(i),colperm1(j)) ...
             - (targ1(i,j) - addcon1);
     
     
     
 end
end
 

if (nopt == 1)
    
    [targ2,vaftarg2,outperm2,rowperm2,colperm2,addcon2,coord2] = uniscaltmac(resprox1,inpermtwo,kblock);
    
end

if (nopt == 0)
    
    outperm2 = inpermtwo;
    [targ2,vaftarg2,rowperm2,colperm2,addcon2,coord2] = linfittmac(resprox1,inpermtwo);
    
end
    


for i = 1:nrow
    for j = 1:ncol
        
          resprox2(rowperm2(i),colperm2(j)) = resprox1(rowperm2(i),colperm2(j)) ...
             - (targ2(i,j) - addcon2);
     
     
     
 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;
addconone = addcon1;
addcontwo= addcon2;

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) - addcontwo);
          
 end
end

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

      [targtwo,vaftwo,rowpermtwo,colpermtwo,addcontwo,coordtwo] = linfittmac(resprox,outpermtwo);

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

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

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

axes = zeros(n,2);

for i = 1:n
    axes(outpermone(i),1) = coordone(i);
    axes(outpermtwo(i),2) = coordtwo(i);
end

