function [find,vaf] = ultrafndtm(proxtm,inpermrow,inpermcol)

% ULTRAFNDTM finds and fits a two-mode ultrametric using 
% iterative projection heuristically on a rectangular proximity
% matrix in the $L_{2}$-norm.
% PROXTM is the input proximity matrix (with a 
% dissimilarity interpretation);
% INPERMROW and INPERMCOL are permutations for the row and column 
% objects that determine the order in which the
% inequality constraints are considered;
% FIND is the found least-squares matrix (with variance-accounted-for
% of VAF) to PROXTM satisfying the ultrametric constraints.

[nrow ncol] = size(proxtm);
work = zeros(nrow*nrow*ncol*ncol,1);
find = proxtm;
cr = 1.0;
iterate = 0;

while (cr >= 1.0e-006)
   cr = 0.0;
   indexll = 0;
   iterate = iterate + 1;
   
  
   
   for ione = 1:(nrow-1)
      for itwo = (ione+1):nrow
         for ithree = 1:(ncol-1)
            for ifour = (ithree+1):ncol
               
               jone = inpermrow(ione);
               jtwo = inpermrow(itwo);
               jthree = inpermcol(ithree);
               jfour = inpermcol(ifour);
            
                         
                                      
            p1 = find(jone,jthree);
            p2 = find(jone,jfour);
            p3 = find(jtwo,jthree);
            p4 = find(jtwo,jfour);
            
            if(iterate <= 100)
            
            find(jone,jthree) = find(jone,jthree) - work(indexll+1);
            find(jone,jfour) = find(jone,jfour) - work(indexll+2);
            find(jtwo,jthree) = find(jtwo,jthree) - work(indexll+3);
            find(jtwo,jfour) = find(jtwo,jfour) - work(indexll+4);
            
            end
         
            
            if(((find(jone,jthree) <= find(jtwo,jthree)) & ...
                  (find(jone,jfour) <= find(jtwo,jthree))) & ...
                  ((find(jone,jthree) <= find(jtwo,jfour)) & ...
                  (find(jone,jfour) <= find(jtwo,jfour)))) 
               
             ave = (find(jtwo,jthree) + find(jtwo,jfour))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = 0.0;
               work(indexll+3) = ave - find(jtwo,jthree); 
               work(indexll+4) = ave - find(jtwo,jfour);
               
               find(jtwo,jthree) = ave;
               find(jtwo,jfour) = ave;
               
               indexll = indexll + 4;
               
                       
               
               
               
               
             elseif(((find(jone,jthree) >= find(jtwo,jthree)) & ...
                  (find(jone,jfour) >= find(jtwo,jthree))) & ...
                  ((find(jone,jthree) >= find(jtwo,jfour)) & ...
                  (find(jone,jfour) >= find(jtwo,jfour)))) 
               
               
               ave = (find(jone,jthree) + find(jone,jfour))/2.0;
               
               
   
               work(indexll+1) = ave - find(jone,jthree);
               work(indexll+2) = ave - find(jone,jfour);
               work(indexll+3) = 0.0; 
               work(indexll+4) = 0.0;
               
               find(jone,jthree) = ave;
               find(jone,jfour) = ave;
               
               indexll = indexll + 4;
               
 
                   
             elseif(((find(jone,jthree) <= find(jone,jfour)) & ...
                  (find(jtwo,jthree) <= find(jone,jfour))) & ...
                  ((find(jone,jthree) <= find(jtwo,jfour)) & ...
                  (find(jtwo,jthree) <= find(jtwo,jfour))))
               
               
               
               
               ave = (find(jone,jfour) + find(jtwo,jfour))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - find(jone,jfour);
               work(indexll+3) = 0.0; 
               work(indexll+4) = ave - find(jtwo,jfour);
               
               find(jone,jfour) = ave;
               find(jtwo,jfour) = ave;
               
               indexll = indexll + 4;
               

                
             elseif(((find(jone,jthree) >= find(jone,jfour)) & ...
                  (find(jtwo,jthree) >= find(jone,jfour))) & ...
                  ((find(jone,jthree) >= find(jtwo,jfour)) & ...
                  (find(jtwo,jthree) >= find(jtwo,jfour))))
               
               ave = (find(jone,jthree) + find(jtwo,jthree))/2.0;
   
               work(indexll+1) = ave - find(jone,jthree);
               work(indexll+2) = 0.0;
               work(indexll+3) = ave - find(jtwo,jthree); 
               work(indexll+4) = 0.0;
               
               find(jone,jthree) = ave;
               find(jtwo,jthree) = ave;
               
               indexll = indexll + 4;
               

                
             elseif(((find(jone,jfour) <= find(jone,jthree)) & ...
                  (find(jtwo, jthree) <= find(jone,jthree))) & ...
                  ((find(jone,jfour) <= find(jtwo,jfour)) & ...
                  (find(jtwo,jthree) <= find(jtwo,jfour))))
               
               
               ave = (find(jone,jthree) + find(jtwo,jfour))/2.0;
   
               work(indexll+1) = ave - find(jone,jthree);
               work(indexll+2) = 0.0;
               work(indexll+3) = 0.0; 
               work(indexll+4) = ave - find(jtwo,jfour);
               
               find(jone,jthree) = ave;
               find(jtwo,jfour) = ave;
               
               indexll = indexll + 4;
               

                
             elseif(((find(jone,jfour) >= find(jone,jthree)) & ...
                  (find(jtwo, jthree) >= find(jone,jthree))) & ...
                  ((find(jone,jfour) >= find(jtwo,jfour)) & ...
                  (find(jtwo,jthree) >= find(jtwo,jfour)))) 
                

               
               
               ave = (find(jone,jfour) + find(jtwo,jthree))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - find(jone,jfour);
               work(indexll+3) = ave - find(jtwo,jthree);
               work(indexll+4) = 0.0;
               
               find(jone,jfour) = ave;
               find(jtwo,jthree) = ave;
               
               indexll = indexll + 4;
                       
            end
            
            cr = cr + abs(p1-find(jone,jthree)) + ...
               abs(p2-find(jone,jfour)) ... 
               + abs(p3-find(jtwo,jthree)) + ...
               abs(p4-find(jtwo,jfour)); 
            
         end
      end
   end
end
end

[fit,vaf] = ultrafittm(proxtm,find);

find = fit;

aveprox = sum(sum(proxtm))/(nrow*(ncol));

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

denom = sum(sum((proxtm- aveprox).^2));

vaf = 1 - (diff/denom);
   
     
                        
               
      
               
   
   
   
   
