function [fit, vaf] = csrobpolish(prox,inperm,targ)



n = size(prox,1);
work = zeros(n*(n-1)*(n-2),1);
fit = prox(inperm,inperm);
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)
         
         p1 = fit(jone,jtwo);
         p2 = fit(jone,jtwo+1);
         
         fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
         fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(indexll+2);
         
         if(((abs(targ(jone,jtwo) - targ(jone,jtwo+1)) < 1.0e-006) ...
               || ((targ(jone,jtwo) < targ(jone,jtwo+1)) && ...
               (fit(jone,jtwo) > fit(jone,jtwo+1)))) || ...
               ((targ(jone,jtwo) > targ(jone,jtwo+1)) && ...
               (fit(jone,jtwo) < fit(jone,jtwo+1))))
            
            ave = (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0;
            work(indexll+1) = ave - fit(jone,jtwo);
            work(indexll+2) = ave - fit(jone,jtwo+1);
            fit(jone,jtwo) = ave;
            fit(jone,jtwo+1) = ave;
            
            cr = cr + abs(fit(jone,jtwo) - p1) + ...
               abs(fit(jone,jtwo+1) - p2);
            
         else
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
            
         end
         
         indexll = indexll + 2;
         
      end
   end
   
          
   for jone = 3:n
      for jtwo = 1:(jone-2)
         
         p1 = fit(jtwo,jone);
         p2 = fit(jtwo+1,jone);
         
         fit(jtwo,jone) = fit(jtwo,jone) - work(indexll+1);
         fit(jtwo+1,jone) = fit(jtwo+1,jone) - work(indexll+2);
         
         if(((abs(targ(jtwo+1,jone) - targ(jtwo,jone)) < 1.0e-006) ...
               || ((targ(jtwo+1,jone) < targ(jtwo,jone)) && ...
               (fit(jtwo+1,jone) > fit(jtwo,jone)))) || ...
               ((targ(jtwo+1,jone) > targ(jtwo,jone)) && ...
               (fit(jtwo+1,jone) < fit(jtwo,jone))))
            
            ave = (fit(jtwo+1,jone) + fit(jtwo,jone))/2.0;
            work(indexll+1) = ave - fit(jtwo,jone);
            work(indexll+2) = ave - fit(jtwo+1,jone);
            fit(jtwo+1,jone) = ave;
            fit(jtwo,jone) = ave;
            
            cr = cr + abs(fit(jtwo,jone) - p1) + ...
               abs(fit(jtwo+1,jone) - p2);
            
         else
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
            
         end
         
         indexll = indexll + 2;
         
      end
   end
   
          
   for jone = 2:(n-1)
   
         
         p1 = fit(jone,n);
         p2 = fit(1,jone);
         
         fit(jone,n) = fit(jone,n) - work(indexll+1);
         fit(1,jone) = fit(1,jone) - work(indexll+2);
         
         if(((abs(targ(jone,n) - targ(1,jone)) < 1.0e-006) ...
               || ((targ(jone,n) < targ(1,jone)) && ...
               (fit(jone,n) > fit(1,jone)))) || ...
               ((targ(jone,n) > targ(1,jone)) && ...
               (fit(jone,n) < fit(1,jone))))
            
            ave = (fit(jone,n) + fit(1,jone))/2.0;
            work(indexll+1) = ave - fit(jone,n);
            work(indexll+2) = ave - fit(1,jone);
            fit(jone,n) = ave;
            fit(1,jone) = ave;
            
            cr = cr + abs(fit(jone,n) - p1) + ...
               abs(fit(1,jone) - p2);
            
         else
            work(indexll+1) = 0.0;
            work(indexll+2) = 0.0;
            
         end
         
         indexll = indexll +2;
         
      end
      
     
      for jone = 1:(n-3)
         for jtwo = (jone+2):(n-1)
            
            p1 = fit(jone,jtwo);
            p2 = fit(jone+1,jtwo);
            p3 = fit(jone,jtwo+1);
            p4 = fit(jone+1,jtwo+1);
            
            if(targ(jone,jtwo+1) <= targ(jone+1,jtwo))
               
               fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
               fit(jone+1,jtwo) = fit(jone+1,jtwo) - work(indexll+2);
               
               if(abs(targ(jone,jtwo+1) - targ(jone+1,jtwo+1)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone,jtwo) + fit(jone+1,jtwo))/2.0;
                  
                  work(indexll+1) = ave - fit(jone,jtwo);
                  work(indexll+2) = ave - fit(jone+1,jtwo);
                  
                  fit(jone,jtwo) = ave;
                  fit(jone+1,jtwo) = ave;
                  
                  cr = cr + abs(fit(jone,jtwo) - p1) + ...
                     abs(fit(jone+1,jtwo) - p2);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            else
               
               fit(jone,jtwo+1) = fit(jone,jtwo+1) - work(indexll+1);
               fit(jone+1,jtwo+1) = fit(jone+1,jtwo+1) - ...
                  work(indexll+2);
               
               if(abs(targ(jone,jtwo) - targ(jone+1,jtwo)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone,jtwo+1) + fit(jone+1,jtwo+1))/2.0;
                  
                  work(indexll+1) = ave - fit(jone,jtwo+1);
                  work(indexll+2) = ave - fit(jone+1,jtwo+1);
                  
                  fit(jone,jtwo+1) = ave;
                  fit(jone+1,jtwo+1) = ave;
                  
                  cr = cr + abs(fit(jone,jtwo+1) - p3) ...
                     + abs(fit(jone+1,jtwo+1) - p4);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            end
         
                
            p1 = fit(jone,jtwo);
            p2 = fit(jone+1,jtwo);
            p3 = fit(jone,jtwo+1);
            p4 = fit(jone+1,jtwo+1);
            
            if(targ(jone,jtwo+1) <= targ(jone+1,jtwo))
               
               fit(jone+1,jtwo) = fit(jone+1,jtwo) - work(indexll+1);
               fit(jone+1,jtwo+1) = fit(jone+1,jtwo+1) - work(indexll+2);
               
               if(abs(targ(jone,jtwo) - targ(jone,jtwo+1)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone+1,jtwo) + fit(jone+1,jtwo+1))/2.0;
                  
                  work(indexll+1) = ave - fit(jone+1,jtwo);
                  work(indexll+2) = ave - fit(jone+1,jtwo+1);
                  
                  fit(jone+1,jtwo) = ave;
                  fit(jone+1,jtwo+1) = ave;
                  
                  cr = cr + abs(fit(jone+1,jtwo) - p2) + ...
                     abs(fit(jone+1,jtwo+1) - p4);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            else
               
               fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
               fit(jone,jtwo+1) = fit(jone,jtwo+1) - ...
                  work(indexll+2);
               
               if(abs(targ(jone+1,jtwo) - targ(jone+1,jtwo+1)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone,jtwo) + fit(jone,jtwo+1))/2.0;
                  
                  work(indexll+1) = ave - fit(jone,jtwo);
                  work(indexll+2) = ave - fit(jone,jtwo+1);
                  
                  fit(jone,jtwo) = ave;
                  fit(jone,jtwo+1) = ave;
                  
                  cr = cr + abs(fit(jone,jtwo) - p1) ...
                     + abs(fit(jone,jtwo+1) - p3);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            end
            
        end
      end
      
                        
      for jone = 2:(n-2)
         
               
            p1 = fit(jone,n);
            p2 = fit(jone+1,n);
            p3 = fit(1,jone);
            p4 = fit(1,jone+1);
            
            if(targ(1,jone) <= targ(jone+1,n))
               
               fit(jone,n) = fit(jone,n) - work(indexll+1);
               fit(jone+1,n) = fit(jone+1,n) - work(indexll+2);
               
               if(abs(targ(1,jone) - targ(1,jone+1)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone,n) + fit(jone+1,n))/2.0;
                  
                  work(indexll+1) = ave - fit(jone,n);
                  work(indexll+2) = ave - fit(jone+1,n);
                  
                  fit(jone,n) = ave;
                  fit(jone+1,n) = ave;
                  
                  cr = cr + abs(fit(jone,n) - p1) + ...
                     abs(fit(jone+1,n) - p2);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            else
               
               fit(1,jone) = fit(1,jone) - work(indexll+1);
               fit(1,jone+1) = fit(1,jone+1) - ...
                  work(indexll+2);
               
               if(abs(targ(jone,n) - targ(jone+1,n)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(1,jone) + fit(1,jone+1))/2.0;
                  
                  work(indexll+1) = ave - fit(1,jone);
                  work(indexll+2) = ave - fit(1,jone+1);
                  
                  fit(1,jone) = ave;
                  fit(1,jone+1) = ave;
                  
                  cr = cr + abs(fit(1,jone) - p3) ...
                     + abs(fit(1,jone+1) - p4);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            end
            
                
            p1 = fit(jone,n);
            p2 = fit(jone+1,n);
            p3 = fit(1,jone);
            p4 = fit(1,jone+1);
            
            if(targ(1,jone) <= targ(jone+1,n))
               
               fit(jone+1,n) = fit(jone+1,n) - work(indexll+1);
               fit(1,jone+1) = fit(1,jone+1) - work(indexll+2);
               
               if(abs(targ(jone,n) - targ(1,jone)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone+1,n) + fit(1,jone+1))/2.0;
                  
                  work(indexll+1) = ave - fit(jone+1,n);
                  work(indexll+2) = ave - fit(1,jone+1);
                  
                  fit(jone+1,n) = ave;
                  fit(1,jone+1) = ave;
                  
                  cr = cr + abs(fit(jone+1,n) - p2) + ...
                     abs(fit(1,jone+1) - p4);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            else
               
               fit(jone,n) = fit(jone,n) - work(indexll+1);
               fit(1,jone) = fit(1,jone) - ...
                  work(indexll+2);
               
               if(abs(targ(jone+1,n) - targ(1,jone+1)) ...
                     <= 1.0e-006)
                  
                  ave = (fit(jone,n) + fit(1,jone))/2.0;
                  
                  work(indexll+1) = ave - fit(jone,n);
                  work(indexll+2) = ave - fit(1,jone);
                  
                  fit(jone,n) = ave;
                  fit(1,jone) = ave;
                  
                  cr = cr + abs(fit(jone,n) - p1) ...
                     + abs(fit(1,jone) - p3);
                  
               else
                  
                  work(indexll+1) = 0.0;
                  work(indexll+2) = 0.0;
                  
               end
               
               indexll = indexll + 2;
               
            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(inperm,inperm) - fit).^2));

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

vaf = 1 - (diff/denom);
