function [fit,vaf] = atreefit(prox,targ)

% ATREEFIT fits a given additive tree using iterative projection to
% a symmetric proximity matrix in the $L_{2}$-norm.
% PROX is the input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation);
% TARG is an matrix of the same size as PROX with entries 
% satisfying the four-point additive tree constraints;
% FIT is the least-squares optimal matrix (with variance-accounted-for
% of VAF) to PROX satisfying the additive tree constraints implicit in TARG.

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

while (cr >= 1.0e-006)
   cr = 0.0;
   indexll = 0;
   
   
   for jone = 1:(n-3)
      for jtwo = (jone+1):(n-2)
         for jthree = (jtwo+1):(n-1)
            for jfour = (jthree+1):n
            
                       
            p1 = fit(jone,jtwo);
            p2 = fit(jone,jthree);
            p3 = fit(jone,jfour);
            p4 = fit(jtwo,jthree);
            p5 = fit(jtwo,jfour);
            p6 = fit(jthree,jfour);
            
            
                       
            fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
            fit(jone,jthree) = fit(jone,jthree) - work(indexll+2);
            fit(jone,jfour) = fit(jone,jfour) - work(indexll+3);
            fit(jtwo,jthree) = fit(jtwo,jthree) - work(indexll+4);
            fit(jtwo,jfour) = fit(jtwo,jfour) - work(indexll+5);
            fit(jthree,jfour) = fit(jthree,jfour) - work(indexll+6);
            
                    
            
            if(((targ(jone,jtwo) + targ(jthree,jfour)) <= ...
                  (targ(jone,jthree) + targ(jtwo,jfour))) & ...
                  ((targ(jone,jtwo) + targ(jthree,jfour)) <=  ...
                  (targ(jone,jfour) + targ(jtwo,jthree)))) 
               
               del = (fit(jone,jthree) + fit(jtwo,jfour) - ...
                  (fit(jone,jfour) + fit(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;
                      
               fit(jone,jthree) = fit(jone,jthree) - del;
               fit(jtwo,jfour) = fit(jtwo,jfour) - del;
               fit(jone,jfour) = fit(jone,jfour) + del;
               fit(jtwo,jthree) = fit(jtwo,jthree) + del;
               
               indexll = indexll + 6;
               
             elseif(((targ(jone,jthree) + targ(jtwo,jfour)) <= ...
                  (targ(jone,jtwo) + targ(jthree,jfour))) & ...
                  ((targ(jone,jthree) + targ(jtwo,jfour)) <=  ...
                  (targ(jone,jfour) + targ(jtwo,jthree)))) 
               
               del = (fit(jone,jtwo) + fit(jthree,jfour) - ...
                  (fit(jone,jfour) + fit(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;
               
               fit(jone,jtwo) = fit(jone,jtwo) - del;
               fit(jthree,jfour) = fit(jthree,jfour) - del;
               fit(jone,jfour) = fit(jone,jfour) + del;
               fit(jtwo,jthree) = fit(jtwo,jthree) + del;
               
               indexll = indexll + 6;
               
            elseif(((targ(jone,jfour) + targ(jtwo,jthree)) <= ...
                  (targ(jone,jtwo) + targ(jthree,jfour))) & ...
                  ((targ(jone,jfour) + targ(jtwo,jthree)) <=  ...
                  (targ(jone,jthree) + targ(jtwo,jfour)))) 
               
               del = (fit(jone,jtwo) + fit(jthree,jfour) - ...
                  (fit(jone,jthree) + fit(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;
               
               fit(jone,jtwo) = fit(jone,jtwo) - del;
               fit(jthree,jfour) = fit(jthree,jfour) - del;
               fit(jone,jthree) = fit(jone,jthree) + del;
               fit(jtwo,jfour) = fit(jtwo,jfour) + del;
               
               indexll = indexll + 6;
            end
            
    
               cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) ... 
                  + abs(p3-fit(jone,jfour)) + abs(p4-fit(jtwo,jthree)) ...
                  + abs(p5-fit(jtwo,jfour)) + abs(p6-fit(jthree,jfour)); 
            end
         end
      end
   end
end

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

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 - fit).^2));

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

vaf = 1 - (diff/denom);

   
     
                        
               
      
               
   
   
   
   
