function [find,vaf,targone,targtwo,outpermone,outpermtwo, ...
    addconone, addcontwo] = bicirac(prox,inperm,kblock)

% BICIRAC finds and fits the sum of two circular
% unidimensional scales using iterative projection to
% a symmetric 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, ...
%     addconone,addcontwo] = bicirac(prox,inperm,kblock)
%
% PROX is the input proximity matrix ($n \times n$ with a zero
% main diagonal and a dissimilarity interpretation);
% INPERM is a given starting permutation of the first $n$ integers;
% FIND is the least-squares optimal matrix (with variance-
% accounted-for of VAF) to PROX and is the sum of the two 
% circular anti-Robinson matrices;
% TARGONE and TARGTWO are based on the two row and column
% object orderings given by the ending permutations OUTPERMONE 
% and OUTPERMTWO. 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.

%tic;
n = size(prox,1);
resprox1 = zeros(n,n);
resprox2 = zeros(n,n);
proxave = zeros(n,n);
resprox = zeros(n,n);

[targ1,vaftarg1,outperm1,addcon1] = unicirac(prox,inperm,kblock);


for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            resprox1(outperm1(i),outperm1(j)) = prox(outperm1(i),outperm1(j)) ...
             - (targ1(i,j) - addcon1);
     else
         resprox1(outperm1(i),outperm1(j)) = 0;
     end
     
 end
end
            
[targ2,vaftarg2,outperm2,addcon2] = unicirac(resprox1,inperm,kblock);

for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            resprox2(outperm2(i),outperm2(j)) = resprox1(outperm2(i),outperm2(j)) ...
             - (targ2(i,j) - addcon2);
     else
         resprox2(outperm2(i),outperm2(j)) = 0;
     end
     
 end
end

find = prox - resprox2;


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.0;
        end
    end
end

diff = sum(sum((prox - find).^2));
denom = sum(sum((prox - proxave).^2));
vaf = 1 - (diff/denom);
            
targone = targ1;
targtwo = targ2;
outpermone = outperm1;
outpermtwo = outperm2;
addconone = addcon1;
addcontwo= addcon2;

vafdiff = 1.0;
iteration = 0;

while ((vafdiff >= 1.0e-006) && (iteration < 10))
    
    vafprev = vaf;
    iteration = iteration + 1;
    
   for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            resprox(outpermtwo(i),outpermtwo(j)) = prox(outpermtwo(i),outpermtwo(j)) ...
             - (targtwo(i,j) - addcontwo);
     else
         resprox(outpermtwo(i),outpermtwo(j)) = 0;
     end
     
 end
end

      [targone,vafone,addconone] = cirfitac_ftarg(resprox,outpermone,targone);
    
  for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            resprox(outpermone(i),outpermone(j)) = prox(outpermone(i),outpermone(j)) ...
             - (targone(i,j) - addconone);
     else
         resprox(outpermone(i),outpermone(j)) = 0;
     end
     
 end
end

    [targtwo,vaftwo,addcontwo] = cirfitac_ftarg(resprox,outpermtwo,targtwo);
   
  for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            find(outpermone(i),outpermone(j)) = targone(i,j) - addconone;
     else
         find(outpermone(i),outpermone(j)) = 0;
     end
     
 end
end

for i = 1:n
    for j = 1:n
        
     if (i ~= j)
            find(outpermtwo(i),outpermtwo(j)) = find(outpermtwo(i),outpermtwo(j)) ...
             + (targtwo(i,j) - addcontwo);
     else
         find(outpermtwo(i),outpermtwo(j)) = 0;
     end
     
 end
end

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

%toc;