function [fit, diff] = cirfit(prox,inperm)

%CIRFIT does a confirmatory fitting of a given order
%  (assumed to reflect a circular ordering around a closed
%  unidimensional structure) using Dykstra's
%  (Kaczmarz's) iterative projection least-squares method.
%  INPERM is the given order; FIT is an $n \times n$ matrix that
%  is fitted to PROX(INPERM,INPERM) with least-squares value DIFF.


n = size(prox,1);
work = zeros(n*(n-1)*(n-2),1);
fit = prox(inperm,inperm);
cr = 1.0;
iterate = 0;
coor = zeros(n,1);

while (cr >= 1.0e-006)
   
  
   cr = 0.0;
   indexll = 0;
   iterate = iterate +1;
  
   
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
         
         p1 = fit(jone,jtwo);
         
         if(iterate <= 10)
            
            fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
            
         end
         
         
         if(fit(jone,jtwo) < 0.0)
            
            work(indexll+1) = -fit(jone,jtwo);
            fit(jone,jtwo) = 0.0;
            
         elseif(fit(jone,jtwo) >= 0.0)
            
            work(indexll+1) = 0;
            
         end
         
         indexll = indexll + 1;
         cr = cr + abs(p1-fit(jone,jtwo));
         
      end
   end
   
   
   for jone = 1:(n-2)
      for jtwo = (jone+2):n
         
         if((jone ~= 1) | (jtwo ~= n))
            
            nhto = jtwo - jone;
            nhaf = n - nhto;
            
            for i = 1:(n-1)
               
               coor(i) = fit(i,i+1);
               
            end
            
            p1 = fit(jone,jtwo);
            coor(n) = fit(1,n);
            
            if(iterate <= 10)
               
               fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+3);
               
               for i = 1:(n-1)
                  if(i < jone)
                     fit(i,i+1) = fit(i,i+1) - work(indexll+2);
                  elseif((i >= jone) & (i < jtwo))
                     fit(i,i+1) = fit(i,i+1) - work(indexll+1);
                  elseif(i >= jtwo)
                     fit(i,i+1) = fit(i,i+1) - work(indexll+2);
                  end
               end
                              
               fit(1,n) = fit(1,n) - work(indexll+2);
            end
            
               
               aa = 0.0;
               bb = 0.0;
               aac = 0.0;
               bbc = 0.0;
               
               for i = 1:(n-1)
                  if(i < jone)
                     bb = bb + fit(i,i+1);
                  elseif((i >= jone) & (i < jtwo))
                     aa = aa + fit(i,i+1);
                  elseif(i >= jtwo)
                     bb = bb + fit(i,i+1);
                  end
               end
               bb = bb + fit(1,n);
               
               if(aa <= bb)
                  del = (aa - fit(jone,jtwo))/(nhto+1);
                  fit(jone,jtwo) = fit(jone,jtwo) + del;
                  for i = 1:nhto
                     fit(jone+i-1,jone+i) = ...
                        fit(jone+i-1,jone+i) - del;
                  end
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = 0.0;
                  work(indexll+3) = del;
                  
                  cr = cr + abs(p1-fit(jone,jtwo));
                  indexll = indexll + 3;
                  
                  
               elseif(aa > bb)
                  
                  del = (bb - fit(jone,jtwo))/(nhaf+1);
                  fit(jone,jtwo) = fit(jone,jtwo) + del;
                  if(jone ~= 1)
                     for i = 1:(jone-1)
                        fit(i,i+1) = fit(i,i+1) - del;
                     end
                  end
                  if(jtwo ~= n)
                     for i = jtwo:(n-1)
                           fit(i,i+1) = fit(i,i+1) - del; 
                     end
                  end
                                 
                        fit(1,n) = fit(1,n) - del;
                        
                        work(indexll+1) = 0.0;
                        work(indexll+2) = -del;
                        work(indexll+3) = del;
                        
                        cr = cr + abs(p1 - fit(jone,jtwo));
                        indexll = indexll +3;
                                       
                  end
             end
         end
      end
   end
   
   targ = fit;
   
work = zeros(n*(n-1)*(n-2),1);
fit = prox(inperm,inperm);
cr = 1.0;
coor = zeros(n,1);

while (cr >= 1.0e-006)
   
   
   cr = 0.0;
   indexll = 0;
     
   
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
         
         p1 = fit(jone,jtwo);
         fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+1);
         
         if(fit(jone,jtwo) < 0.0)
            
            work(indexll+1) = -fit(jone,jtwo);
            fit(jone,jtwo) = 0.0;
            
         elseif(fit(jone,jtwo) >= 0.0)
            
            work(indexll+1) = 0.0;
            
         end
         
         indexll = indexll + 1;
         cr = cr + abs(p1-fit(jone,jtwo));
         
      end
   end
   
   
   for jone = 1:(n-2)
      for jtwo = (jone+2):n
         
         if((jone ~= 1) | (jtwo ~= n))
            
            nhto = jtwo - jone;
            nhaf = n - nhto;
            
            for i = 1:(n-1)
               
               coor(i) = fit(i,i+1);
               
            end
            
            p1 = fit(jone,jtwo);
            coor(n) = fit(1,n);
            
                  
               fit(jone,jtwo) = fit(jone,jtwo) - work(indexll+3);
               
               for i = 1:(n-1)
                  if(i < jone)
                     fit(i,i+1) = fit(i,i+1) - work(indexll+2);
                  elseif((i >= jone) & (i < jtwo))
                     fit(i,i+1) = fit(i,i+1) - work(indexll+1);
                  elseif(i >= jtwo)
                     fit(i,i+1) = fit(i,i+1) - work(indexll+2);
                  end
               end
               
                  fit(1,n) = fit(1,n) - work(indexll+2);
               
   
            
               
               aa = 0.0;
               bb = 0.0;
               aac = 0.0;
               bbc = 0.0;
               
               for i = 1:(n-1)
                  if(i < jone)
                     bbc = bbc + targ(i,i+1);
                  elseif((i >= jone) & (i < jtwo))
                     aac = aac + targ(i,i+1);
                  elseif(i >= jtwo)
                     bbc = bbc + targ(i,i+1);
                  end
               end
               
               bbc = bbc + targ(1,n);
               
               
               for i = 1:(n-1)
                  if(i < jone)
                     bb = bb + fit(i,i+1);
                  elseif((i >= jone) & (i < jtwo))
                     aa = aa + fit(i,i+1);
                  elseif(i >= jtwo)
                     bb = bb + fit(i,i+1);
                  end
               end
               bb = bb + fit(1,n);
               
               if(aac <= bbc)
                  del = (aa - fit(jone,jtwo))/(nhto+1);
                  fit(jone,jtwo) = fit(jone,jtwo) + del;
                  for i = 1:nhto
                     fit(jone+i-1,jone+i) = ...
                        fit(jone+i-1,jone+i) - del;
                  end
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = 0.0;
                  work(indexll+3) = del;
                  
                  cr = cr + abs(p1-fit(jone,jtwo));
                  indexll = indexll + 3;
                  
                  
               elseif(aac > bbc)
                  del = (bb - fit(jone,jtwo))/(nhaf+1);
                  fit(jone,jtwo) = fit(jone,jtwo) + del;
                  if(jone ~= 1)
                     for i = 1:(jone-1)
                        fit(i,i+1) = fit(i,i+1) - del;
                     end
                  end
                  if(jtwo ~= n)
                     for i = jtwo:(n-1)
                        fit(i,i+1) = fit(i,i+1) - del; 
                     end
                  end
                                    
                        fit(1,n) = fit(1,n) - del;
                        
                        work(indexll+1) = 0.0;
                        work(indexll+2) = -del;
                        work(indexll+3) = del;
                        
                        cr = cr + abs(p1 - fit(jone,jtwo));
                        indexll = indexll +3;
                      
                     end
             end
         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)).^2));

vaf = 1 - (diff/denom);
  
diff = (.5)*diff;
   
   
   
 
   
         
                        
                        
                     
                     
                     
                     
