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

% 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.
%
% syntax: [find,vaf,targone,targtwo,outpermone,outpermtwo, ...
%        rowpermone,colpermone,rowpermtwo,colpermtwo,addconone,...
%        addcontwo,coordone,coordtwo,axes] = ...
%    biscaltmac(proxtm,inpermone,inpermtwo,kblock,nopt)
%
% PROXTM is the input two-mode proximity matrix ($nrow \times ncol$
% with 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 of 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.

%tic;
[nrow ncol] = size(proxtm);
n = nrow + ncol;
resprox1 = zeros(nrow,ncol);
resprox2 = zeros(nrow,ncol);
proxave = zeros(nrow,ncol);
resprox = zeros(nrow,ncol);
%rowdum = zeros(nrow,1);
%coldum = zeros(ncol,1);
%dumone = zeros(1,1);
invar = zeros(2,1);


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

      [targ1,vaftarg1,outperm1,rowperm1,colperm1,addcon1,coord1] =  ms_uniscaltmac(proxtm,inpermone,kblock,addcontol,maxiteration);
      
end

if (nopt == 0)
    
    outperm1 = inpermone;
%    [targ1,vaftarg1,rowperm1,colperm1,addcon1,coord1] = ms_linfittmac(proxtm,inpermone,rowdum,coldum);
    
    [targ1,coord1,outvar] = ms_linfittmac(proxtm,inpermone,invar,addcontol,maxiteration);
    
    vaftarg1 = outvar(1);
    addcon1 = outvar(2);


   rowperm1 = zeros(nrow,1);
   colperm1 = zeros(ncol,1);
   colidx = 0;
   rowidx = 0;
   
   for i = 1:n
      
      if (inpermone(i) <= nrow)
         
         rowidx = rowidx + 1;
        rowperm1(rowidx) = inpermone(i);
      end
         
      if (inpermone(i) > nrow)
         
         colidx = colidx + 1;
         colperm1(colidx) = inpermone(i) - nrow;
         
      end
   end
    
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] = ms_uniscaltmac(resprox1,inpermtwo,kblock);
    
    [targ2,vaftarg2,outperm2,rowperm2,colperm2,addcon2,coord2] =  ms_uniscaltmac(resprox1,inpermtwo,kblock,addcontol,maxiteration);
    
 
    
end

if (nopt == 0)
    
    outperm2 = inpermtwo;
%    [targ2,vaftarg2,rowperm2,colperm2,addcon2,coord2] = ms_linfittmac(resprox1,inpermtwo,rowdum,coldum,dumone);
    
    
    [targ2,coord2,outvar] = ms_linfittmac(resprox1,inpermtwo,invar,addcontol,maxiteration);
    
    vaftarg2 = outvar(1);
    addcon2 = outvar(2);


   rowperm2 = zeros(nrow,1);
   colperm2 = zeros(ncol,1);
   colidx = 0;
   rowidx = 0;
   
   for i = 1:n
      
      if (inpermtwo(i) <= nrow)
         
         rowidx = rowidx + 1;
        rowperm2(rowidx) = inpermtwo(i);
      end
         
      if (inpermtwo(i) > nrow)
         
         colidx = colidx + 1;
         colperm2(colidx) = inpermtwo(i) - nrow;
         
      end
   end
    
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] = ms_linfittmac(resprox,outpermone,rowdum,coldum,dumone);
    
          [targone,coordone,outvar] = ms_linfittmac(resprox,outpermone,invar,addcontol,maxiteration);
          
          vafone = outvar(1);
          addconone = outvar(2);


   rowpermone= zeros(nrow,1);
   colpermone = zeros(ncol,1);
   colidx = 0;
   rowidx = 0;
   
   for i = 1:n
      
      if (outpermone(i) <= nrow)
         
         rowidx = rowidx + 1;
        rowpermone(rowidx) = outpermone(i);
      end
         
      if (outpermone(i) > nrow)
         
         colidx = colidx + 1;
         colpermone(colidx) = outpermone(i) - nrow;
         
      end
   end
          
  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] = ms_linfittmac(resprox,outpermtwo,rowdum,coldum,dumone);

      [targtwo,coordtwo,outvar] = ms_linfittmac(resprox,outpermtwo,invar,addcontol,maxiteration);
      
      
      vaftwo = outvar(1);
      addcontwo = outvar(2);


   rowpermtwo = zeros(nrow,1);
   colpermtwo = zeros(ncol,1);
   colidx = 0;
   rowidx = 0;
   
   for i = 1:n
      
      if (outpermtwo(i) <= nrow)
         
         rowidx = rowidx + 1;
        rowpermtwo(rowidx) = outpermtwo(i);
      end
         
      if (outpermtwo(i) > nrow)
         
         colidx = colidx + 1;
         colpermtwo(colidx) = outpermtwo(i) - nrow;
         
      end
   end
      
  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

%toc;

