function [find,vaf] = ultrafnd_missing(prox,inperm,proxmiss)

% ULTRAFND_MISSING finds and fits an ultrametric using iterative projection
% heuristically on a symmetric proximity matrix in the $L_{2}$-norm.
%
% syntax: [find,vaf] = ultrafnd_missing(prox,inperm,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.
% 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 ultrametric constraints. 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-2),1);
find = prox;
cr = 1.0;
iterate = 0;
proxave = zeros(n,n);

while (cr >= 1.0e-006)
   cr = 0.0;
   indexll = 0;
   iterate = iterate + 1;
   
   
   for ione = 1:(n-2)
      for itwo = (ione+1):(n-1)
        for ithree = (itwo+1):n
            
            indxord(1) = inperm(ione);
            indxord(2) = inperm(itwo);
            indxord(3) = inperm(ithree);
            
            for i = 1:2
               for j = (i+1):3
                  
                  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);
            
            if(((proxmiss(jone,jtwo) == 1) & (proxmiss(jone,jthree) == 1)) ...
                     & (proxmiss(jtwo,jthree) == 1))
                          
            p1 = find(jone,jtwo);
            p2 = find(jone,jthree);
            p3 = find(jtwo,jthree);
            
            if(iterate <= 100)
            
            find(jone,jtwo) = find(jone,jtwo) - work(indexll+1);
            find(jone,jthree) = find(jone,jthree) - work(indexll+2);
            find(jtwo,jthree) = find(jtwo,jthree) - work(indexll+3);
            
            end
         
            
            if((find(jone,jtwo) <= find(jone,jthree)) && ...
                  (find(jone,jtwo) <= find(jtwo,jthree))) 
               
               ave = (find(jone,jthree) + find(jtwo,jthree))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - find(jone,jthree);
               work(indexll+3) = ave - find(jtwo,jthree); 
               
               find(jone,jthree) = ave;
               find(jtwo,jthree) = ave;
               
               indexll = indexll + 3;
                       
             elseif((find(jone,jthree) <= find(jone,jtwo)) && ...
                  (find(jone,jthree) <= find(jtwo,jthree))) 
               
               ave = (find(jone,jtwo) + find(jtwo,jthree))/2.0;
   
               work(indexll+1) = ave - find(jone,jtwo);
               work(indexll+2) = 0.0;
               work(indexll+3) = ave - find(jtwo,jthree); 
               
               find(jone,jtwo) = ave;
               find(jtwo,jthree) = ave;
               
               indexll = indexll + 3;
   
             elseif((find(jtwo,jthree) <= find(jone,jthree)) && ...
                  (find(jtwo,jthree) <= find(jone,jtwo))) 
               
               ave = (find(jone,jthree) + find(jone,jtwo))/2.0;
   
               work(indexll+1) = ave - find(jone,jtwo);
               work(indexll+2) = ave - find(jone,jthree);
               work(indexll+3) = 0.0; 
               
               find(jone,jtwo) = ave;
               find(jone,jthree) = ave;
               
               indexll = indexll + 3;
            end
            
               cr = cr + abs(p1-find(jone,jtwo)) + abs(p2-find(jone,jthree)) ... 
               + abs(p3-find(jtwo,jthree)); 
           
            end
            
         end
      end
   end
end

  npres = 0;
    
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
                  
         find(jtwo,jone) = find(jone,jtwo);
         
         if(proxmiss(jone,jtwo) == 1)
             
             npres = npres + 1;
         end
         
      end
   end
   
   npres = 2*npres;
   
 [fit, vaf] = ultrafit_missing(prox,find,proxmiss);

find = fit;

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

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

vaf = 1 - (diff/denom);
      


   
     
                        
               
      
               
   
   
   
   
