function [fit,vaf] = ultrafit_missing(prox,targ,proxmiss)

% ULTRAFIT_MISSING fits a given ultrametric using iterative projection to
% a symmetric proximity matrix in the $L_{2}$-norm.
%
% syntax: [fit,vaf] = ultrafit_missing(prox,targ,proxmiss)
%
% PROX is the input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation); also, missing entries in the input
% proximity matrix PROX are given values of zero.
% TARG is an ultrametric matrix of the same size as PROX;
% FIT is the least-squares optimal matrix (with 
% variance-accounted-for of VAF) to PROX satisfying the ultrametric
% constraints implicit in TARG. PROXMISS is the same size as PROX (with main
% diagonal entries all zero); an off-diagonal entry of 1.0 denotes an
% entry in PROX that is present and 0.0 if it is absent.

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

while (cr >= 1.0e-006)
   cr = 0.0;
   indexll = 0;
   
   for jone = 1:(n-2)
      for jtwo = (jone+1):(n-1)
         for jthree = (jtwo+1):n
             
             if(((proxmiss(jone,jtwo) == 1) & (proxmiss(jone,jthree) == 1)) ...
                     & (proxmiss(jtwo,jthree) == 1))
            
            p1 = fit(jone,jtwo);
            p2 = fit(jone,jthree);
            p3 = fit(jtwo,jthree);
            
            fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
            fit(jone,jthree) = fit(jone,jthree) - work(indexll+2);
            fit(jtwo,jthree) = fit(jtwo,jthree) - work(indexll+3);
            
            if((targ(jone,jtwo) <= targ(jone,jthree)) && ...
                  (targ(jone,jtwo) <= targ(jtwo,jthree))) 
               
               ave = (fit(jone,jthree) + fit(jtwo,jthree))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - fit(jone,jthree);
               work(indexll+3) = ave - fit(jtwo,jthree); 
               
               fit(jone,jthree) = ave;
               fit(jtwo,jthree) = ave;
               
               indexll = indexll + 3;
                       
             elseif((targ(jone,jthree) <= targ(jone,jtwo)) && ...
                  (targ(jone,jthree) <= targ(jtwo,jthree))) 
               
               ave = (fit(jone,jtwo) + fit(jtwo,jthree))/2.0;
   
               work(indexll+1) = ave - fit(jone,jtwo);
               work(indexll+2) = 0.0;
               work(indexll+3) = ave - fit(jtwo,jthree); 
               
               fit(jone,jtwo) = ave;
               fit(jtwo,jthree) = ave;
               
               indexll = indexll + 3;
   
             elseif((targ(jtwo,jthree) <= targ(jone,jthree)) && ...
                  (targ(jtwo,jthree) <= targ(jone,jtwo))) 
               
               ave = (fit(jone,jthree) + fit(jone,jtwo))/2.0;
   
               work(indexll+1) = ave - fit(jone,jtwo);
               work(indexll+2) = ave - fit(jone,jthree);
               work(indexll+3) = 0.0; 
               
               fit(jone,jtwo) = ave;
               fit(jone,jthree) = ave;
               
               indexll = indexll + 3;
            end
            
               cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) ... 
               + abs(p3-fit(jtwo,jthree));
             end
            
         end
      end
   end
   
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
         for jthree = 1:(n-1)
            for jfour = (jthree+1):n
                
                if((proxmiss(jone,jtwo) == 1) & (proxmiss(jthree,jfour) == 1))
               
               if((jone ~= jthree) || (jtwo ~= jfour))
                  
                  p1 = fit(jone,jtwo);
                  p2 = fit(jthree,jfour);
                  fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
                  fit(jthree,jfour) = fit(jthree,jfour) - work(indexll+2);
                                    
                  if(abs(targ(jone,jtwo) - targ(jthree,jfour)) <= ...
                        1.0e-006)
                     
                     ave = (fit(jone,jtwo) + fit(jthree,jfour))/2.0;
                     work(indexll+1) = ave - fit(jone,jtwo);
                     work(indexll+2) = ave - fit(jthree,jfour);
                     fit(jone,jtwo) = ave;
                     fit(jthree,jfour) = ave;
                     
                     indexll = indexll + 2;
                     
                  elseif((abs(targ(jone,jtwo) - targ(jthree,jfour)) ...
                        > 1.0e-006) && (targ(jone,jtwo) < ...
                        targ(jthree,jfour)))
                     
                     if(fit(jone,jtwo) < fit(jthree,jfour))
                        
                        work(indexll+1) = 0;
                        work(indexll+2) = 0;
                        
                        indexll = indexll + 2;
                        
                     elseif(fit(jone,jtwo) >= fit(jthree,jfour))
                        
                        ave = (fit(jone,jtwo) + fit(jthree,jfour))/2.0;
                        work(indexll+1) = ave - fit(jone,jtwo);
                        work(indexll+2) = ave - fit(jthree,jfour);
                        fit(jone,jtwo) = ave;
                        fit(jthree,jfour) = ave;
                        
                        indexll = indexll + 2;
                        
                     end
                     
                  elseif((abs(targ(jone,jtwo) - targ(jthree,jfour)) ...
                        > 1.0e-006) && (targ(jone,jtwo) > ...
                        targ(jthree,jfour)))
                     
                     if(fit(jone,jtwo) > fit(jthree,jfour))
                        
                        work(indexll+1) = 0;
                        work(indexll+2) = 0;
                        
                        indexll = indexll + 2;
                        
                     elseif(fit(jone,jtwo) <= fit(jthree,jfour))
                        
                        ave = (fit(jone,jtwo) + fit(jthree,jfour))/2.0;
                        work(indexll+1) = ave - fit(jone,jtwo);
                        work(indexll+2) = ave - fit(jthree,jfour);
                        fit(jone,jtwo) = ave;
                        fit(jthree,jfour) = ave;
                        
                        indexll = indexll + 2;
                        
                     end
                  end
                             
               cr = cr + abs(p1-fit(jone,jtwo)) + ...
                  abs(p2-fit(jthree,jfour));
               end
                end
            
               
            end
         end
      end
   end
end

   
      
                  
    npres = 0;
    
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
                  
         fit(jtwo,jone) = fit(jone,jtwo);
         
         if(proxmiss(jone,jtwo) == 1)
             
             npres = npres + 1;
         end
         
      end
   end
   
   npres = 2*npres;
   
  aveprox = sum(sum(prox))/(npres);

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

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

vaf = 1 - (diff/denom);
      
               
   
   
   
   
