function [find,vaf] = atreefnd(prox,inperm)

% ATREEFIND finds and fits an additive tree using iterative projection
% heuristically on a symmetric proximity matrix in the $L_{2}$-norm.
% 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 satisfying the additive tree constraints.


n = size(prox,1);
work = zeros(n*(n-1)*(n-2)*(n-3),1);
find = prox;
cr = 1.0;
iterate = 0;

while (cr >= 1.0e-006)
   cr = 0.0;
   indexll = 0;
   iterate = iterate + 1;
   
   
   for ione = 1:(n-3)
      for itwo = (ione+1):(n-2)
         for ithree = (itwo+1):(n-1)
            for ifour = (ithree+1):n
            
            indxord(1) = inperm(ione);
            indxord(2) = inperm(itwo);
            indxord(3) = inperm(ithree);
            indxord(4) = inperm(ifour);
            
            for i = 1:3
               for j = (i+1):4
                  
                  if(indxord(i) > indxord(j))
                     
                     temp = indxord(i);
                     indxord(i) = indxord(j);
                     indxord(j) = temp;
                     
                  end
               end
            end
            
               
            jone = indxord(1);
            jtwo = indxord(2);
            jthree = indxord(3);
            jfour = indxord(4);
            
            p1 = find(jone,jtwo);
            p2 = find(jone,jthree);
            p3 = find(jone,jfour);
            p4 = find(jtwo,jthree);
            p5 = find(jtwo,jfour);
            p6 = find(jthree,jfour);
            
            
            if(iterate <= 100)
            
            find(jone,jtwo) = find(jone,jtwo) - work(indexll+1);
            find(jone,jthree) = find(jone,jthree) - work(indexll+2);
            find(jone,jfour) = find(jone,jfour) - work(indexll+3);
            find(jtwo,jthree) = find(jtwo,jthree) - work(indexll+4);
            find(jtwo,jfour) = find(jtwo,jfour) - work(indexll+5);
            find(jthree,jfour) = find(jthree,jfour) - work(indexll+6);
            
            end
         
            
            if(((find(jone,jtwo) + find(jthree,jfour)) <= ...
                  (find(jone,jthree) + find(jtwo,jfour))) & ...
                  ((find(jone,jtwo) + find(jthree,jfour)) <=  ...
                  (find(jone,jfour) + find(jtwo,jthree)))) 
               
               del = (find(jone,jthree) + find(jtwo,jfour) - ...
                  (find(jone,jfour) + find(jtwo,jthree)))/4.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = -del;
               work(indexll+3) =  del;
               work(indexll+4) = del;
               work(indexll+5) = -del;
               work(indexll+6) = 0.0;
                      
               find(jone,jthree) = find(jone,jthree) - del;
               find(jtwo,jfour) = find(jtwo,jfour) - del;
               find(jone,jfour) = find(jone,jfour) + del;
               find(jtwo,jthree) = find(jtwo,jthree) + del;
               
               indexll = indexll + 6;
               
             elseif(((find(jone,jthree) + find(jtwo,jfour)) <= ...
                  (find(jone,jtwo) + find(jthree,jfour))) & ...
                  ((find(jone,jthree) + find(jtwo,jfour)) <=  ...
                  (find(jone,jfour) + find(jtwo,jthree)))) 
               
               del = (find(jone,jtwo) + find(jthree,jfour) - ...
                  (find(jone,jfour) + find(jtwo,jthree)))/4.0;
   
               work(indexll+1) = -del;
               work(indexll+2) = 0.0;
               work(indexll+3) =  del;
               work(indexll+4) = del;
               work(indexll+5) = 0.0;
               work(indexll+6) = -del;
               
               find(jone,jtwo) = find(jone,jtwo) - del;
               find(jthree,jfour) = find(jthree,jfour) - del;
               find(jone,jfour) = find(jone,jfour) + del;
               find(jtwo,jthree) = find(jtwo,jthree) + del;
               
               indexll = indexll + 6;
               
            elseif(((find(jone,jfour) + find(jtwo,jthree)) <= ...
                  (find(jone,jtwo) + find(jthree,jfour))) & ...
                  ((find(jone,jfour) + find(jtwo,jthree)) <=  ...
                  (find(jone,jthree) + find(jtwo,jfour)))) 
               
               del = (find(jone,jtwo) + find(jthree,jfour) - ...
                  (find(jone,jthree) + find(jtwo,jfour)))/4.0;
   
               work(indexll+1) = -del;
               work(indexll+2) =  del;
               work(indexll+3) =  0.0;
               work(indexll+4) =  0.0;
               work(indexll+5) =  del;
               work(indexll+6) = -del;
               
               find(jone,jtwo) = find(jone,jtwo) - del;
               find(jthree,jfour) = find(jthree,jfour) - del;
               find(jone,jthree) = find(jone,jthree) + del;
               find(jtwo,jfour) = find(jtwo,jfour) + del;
               
               indexll = indexll + 6;
            end
            
    
               cr = cr + abs(p1-find(jone,jtwo)) + abs(p2-find(jone,jthree)) ... 
                  + abs(p3-find(jone,jfour)) + abs(p4-find(jtwo,jthree)) ...
                  + abs(p5-find(jtwo,jfour)) + abs(p6-find(jthree,jfour)); 
            end
         end
      end
   end
end

for jone = 1:(n-1)
   for jtwo = (jone+1):n
      
      find(jtwo,jone) = find(jone,jtwo);
      
   end
end


[fit,vaf] = atreefit(prox,find);

find = 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;
      end
   end
end

diff = sum(sum((prox - find).^2));

denom = sum(sum((prox - proxave).^2));

vaf = 1 - (diff/denom);
   
     
                        
               
      
               
   
   
   
   
