function [fit,vaf] = ultrafittm(proxtm,targ)


% ULTRAFITTM fits a given (two-mode) ultrametric using iterative projection to
% a two-mode (rectangular) proximity matrix in the $L_{2}$-norm.
% PROXTM is the input proximity matrix (with a dissimilarity interpretation);
% TARG is an ultrametric matrix of the same size as PROXTM;
% FIT is the least-squares optimal matrix (with variance-accounted-for
% of VAF) to PROXTM satisfying the ultrametric constraints implicit in TARG.

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

while (cr >= 1.0e-006)
   
   
   cr = 0.0;
   indexll = 0;
   iterate = iterate + 1;
    
   
   for jone = 1:(nrow-1)
      for jtwo = (jone+1):nrow
         for jthree = 1:(ncol-1)
            for jfour = (jthree+1):ncol
               
                                        
                                      
            p1 = fit(jone,jthree);
            p2 = fit(jone,jfour);
            p3 = fit(jtwo,jthree);
            p4 = fit(jtwo,jfour);
            
        
            
            fit(jone,jthree) = fit(jone,jthree) - work(indexll+1);
            fit(jone,jfour) = fit(jone,jfour) - work(indexll+2);
            fit(jtwo,jthree) = fit(jtwo,jthree) - work(indexll+3);
            fit(jtwo,jfour) = fit(jtwo,jfour) - work(indexll+4);
            
           
         
            
            if(((targ(jone,jthree) <= targ(jtwo,jthree)) & ...
                  (targ(jone,jfour) <= targ(jtwo,jthree))) & ...
                  ((targ(jone,jthree) <= targ(jtwo,jfour)) & ...
                  (targ(jone,jfour) <= targ(jtwo,jfour)))) 
               
             ave = (fit(jtwo,jthree) + fit(jtwo,jfour))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = 0.0;
               work(indexll+3) = ave - fit(jtwo,jthree); 
               work(indexll+4) = ave - fit(jtwo,jfour);
               
               fit(jtwo,jthree) = ave;
               fit(jtwo,jfour) = ave;
               
               indexll = indexll + 4;
               
                       
               
               
               
               
             elseif(((targ(jone,jthree) >= targ(jtwo,jthree)) & ...
                  (targ(jone,jfour) >= targ(jtwo,jthree))) & ...
                  ((targ(jone,jthree) >= targ(jtwo,jfour)) & ...
                  (targ(jone,jfour) >= targ(jtwo,jfour)))) 
               
               
               ave = (fit(jone,jthree) + fit(jone,jfour))/2.0;
               
               
   
               work(indexll+1) = ave - fit(jone,jthree);
               work(indexll+2) = ave - fit(jone,jfour);
               work(indexll+3) = 0.0; 
               work(indexll+4) = 0.0;
               
               fit(jone,jthree) = ave;
               fit(jone,jfour) = ave;
               
               indexll = indexll + 4;
               
 
                   
             elseif(((targ(jone,jthree) <= targ(jone,jfour)) & ...
                  (targ(jtwo,jthree) <= targ(jone,jfour))) & ...
                  ((targ(jone,jthree) <= targ(jtwo,jfour)) & ...
                  (targ(jtwo,jthree) <= targ(jtwo,jfour))))
               
               
               
               
               ave = (fit(jone,jfour) + fit(jtwo,jfour))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - fit(jone,jfour);
               work(indexll+3) = 0.0; 
               work(indexll+4) = ave - fit(jtwo,jfour);
               
               fit(jone,jfour) = ave;
               fit(jtwo,jfour) = ave;
               
               indexll = indexll + 4;
               

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

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

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

               
               
               ave = (fit(jone,jfour) + fit(jtwo,jthree))/2.0;
   
               work(indexll+1) = 0.0;
               work(indexll+2) = ave - fit(jone,jfour);
               work(indexll+3) = ave - fit(jtwo,jthree);
               work(indexll+4) = 0.0;
               
               fit(jone,jfour) = ave;
               fit(jtwo,jthree) = ave;
               
               indexll = indexll + 4;
               
            else
               
                       

              
               
               
               
            end
            
            cr = cr + abs(p1-fit(jone,jthree)) + ...
               abs(p2-fit(jone,jfour)) ... 
               + abs(p3-fit(jtwo,jthree)) + ...
               abs(p4-fit(jtwo,jfour)); 
            
         end
      end
   end
end






   for jone = 1:nrow
      for jtwo = 1:ncol
         for jthree = jone:nrow
            for jfour = jtwo:ncol
               
               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





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


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

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

vaf = 1 - (diff/denom);
   
     
                        
               
      
               
   
   
   
   
