function [find,vaf,outperm,targone,targtwo,lengthsone] = ... 
    cent_ultrafnd_confnd(prox,inperm)

% CENT_ULTRAFND_CONFND finds and fits an additive tree by first fitting
% a centroid metric and secondly an ultrametric to the residual
% matrix where the latter is displayed by a constaining object order that
% is also identified in the process.
%
% syntax: [find,vaf,outperm,targone,targtwo,lengthsone] = ...
%      cent_ultrafnd_confnd(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 in identifying the ultrametric;
% FIND is the found least-squares matrix (with variance-accounted-for
% of VAF) to PROX satisfying the additive tree constraints. TARGTWO is 
% the ultrametric component of the decomposition; TARGONE is the centroid
% metric component defined by the lengths in LENGTHSONE; OUTPERM is the
% identified constraining object order used to display the ultrametric
% component.


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

[ctfnd,ctvaf,lengths] = centfit(prox);
resprox = prox - ctfnd;

[conperm,rawindex,allperms,index] = order(resprox,targlin(n),inperm,2);

proxpermut = prox(conperm,conperm);

outperm = conperm;

[ctfnd,ctvaf,lengths] = centfit(proxpermut);

resprox = proxpermut - ctfnd;

[fit,vaf,vafarob,arobprox,vafultra] = ultrafnd_confit(resprox,inperm,1:n);

find = ctfnd + fit;

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((proxpermut - find).^2));
denom = sum(sum((proxpermut - proxave).^2));
vaf = 1 - (diff/denom);
            
targone = ctfnd;
targtwo = fit;

vafdiff = 1.0;
iteration = 0;

while ((vafdiff >= 1.0e-006) && (iteration < 100 ))
    
    vafprev = vaf;
    iteration = iteration + 1;
    
    resprox = proxpermut - targtwo;
    [targone,vafone,lengthsone] = centfit(resprox);
    
    resprox = proxpermut - targone;
    [targtwo,vaftwo,vafarob,arobprox,vafultra] = ultrafnd_confit(resprox,inperm,1:n);
    
       
    find = targone + targtwo;
    
    diff = sum(sum((proxpermut - find).^2));
    denom = sum(sum((proxpermut - proxave).^2));
    vaf = 1 - (diff/denom);
    
    vafdiff = abs(vaf - vafprev);
end