function [find,vaf,targone,targtwo] = biatreefnd(prox,inperm)

% BIATREEFND finds and fits the sum 
% of two additive trees using iterative projection
% heuristically on a symmetric proximity matrix in the $L_{2}$-norm.
%
% syntax: [find,vaf,targone,targtwo] = biatreefnd(prox,inperm)
%
% PROX is the input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation);
% INPERM is a permutation that determines the order in which the
% inequality constraints are considered;
% FIND is the found least-squares matrix (with variance-accounted-for
% of VAF) to PROX and is the sum of 
% the two additive tree matrices TARGONE and TARGTWO.


n = size(prox,1);
proxave = zeros(n,n);

[targ1,vaftarg1] = atreefnd(prox,inperm);

resprox = prox - targ1;

[targ2,vaftarg2] = atreefnd(resprox,inperm);

find = targ1 + targ2;

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;

vafdiff = 1.0;
iteration = 0;

while ((vafdiff >= 1.0e-006) && (iteration < 100))
    
    vafprev = vaf;
    iteration = iteration + 1;
    
    resprox = prox - targtwo;
    [targone,vafone] = atreefnd(resprox,inperm);  
    
    
    resprox = prox - targone;
    [targtwo,vaftwo] = atreefnd(resprox,inperm);
    
    find = targone + targtwo;
    
    diff = sum(sum((prox - find).^2));
    denom = sum(sum((prox - proxave).^2));
    vaf = 1 - (diff/denom);
    
    vafdiff = abs(vaf - vafprev);
end