function [find,vaf,targone,targtwo,outpermone,outpermtwo] = bicirarobfnd(prox,inperm,kblock)

% BICIRAROBFND finds and fits the sum of two circular anti-Robinson 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.
% 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 based on the two row and column
% object orderings given by the ending permutations OUTPERMONE 
% and OUTPERMTWO. 

n = size(prox,1);

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

resprox1(outperm1,outperm1) = prox(outperm1,outperm1) - targ1;

[targ2,vaftarg2,outperm2,addcon2] = unicirac(resprox1,inperm,kblock);

resprox2(outperm2,outperm2) = resprox1(outperm2,outperm2) - targ2;

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;

vafdiff = 1.0;

while (vafdiff >= 1.0e-006)
    
    vafprev = vaf;
    
   
    resprox(outpermtwo,outpermtwo) = prox(outpermtwo,outpermtwo) - targtwo;
    [targone,vafone] = cirarobfit(resprox,outpermone,targone);
    
    resprox(outpermone,outpermone) = prox(outpermone,outpermone) - targone;
    [targtwo,vaftwo] = cirarobfit(resprox,outpermtwo,targtwo);
    
    find(outpermone,outpermone) = targone;
    find(outpermtwo,outpermtwo) = find(outpermtwo,outpermtwo) + targtwo;
    
    diff = sum(sum((prox - find).^2));
    denom = sum(sum((prox - proxave).^2));
    vaf = 1 - (diff/denom);
    
    vafdiff = abs(vaf - vafprev);
end