function [fit, vaf, addcon] = ms_cirfitac(prox,inperm,dumone)

%  CIRFITAC does a confirmatory fitting (including
%  the estimation of an additive constant) for a given order
%  (assumed to reflect a circular ordering around a closed
%  unidimensional structure) using the Dykstra--Kaczmarz
%  iterative projection least-squares method.
%
%  syntax: [fit, vaf, addcon] = cirfitac(prox,inperm)
%
%  INPERM is the given order; FIT is an $n \times n$ matrix that
%  is fitted to PROX(INPERM,INPERM) with variance-accounted-for of 
%  VAF; ADDCON is the estimated additive constant.


%tic;

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);
proxave = zeros(n,n);

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;
            
         end
            
         if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     fit(i,i+1) = fit(i,i+1) - work(indexll+1);
                  end
                  if (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;
              
               for i = 1:(n-1)
                  if (i < jone)
                     bb = bb + fit(i,i+1);
                  end
                  if ((i >= jone) & (i < jtwo))
                     aa = aa + fit(i,i+1);
                  end
                  if (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;
               end
                  
                  
         if (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;
         end
            
         if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     fit(i,i+1) = fit(i,i+1) - work(indexll+1);
                  end
                  if (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;
               bbc = 0.0;
               aac = 0.0;
               
               for i = 1:(n-1)
                  if (i < jone)
                     bbc = bbc + targ(i,i+1);
                  end
                  if ((i >= jone) & (i < jtwo))
                     aac = aac + targ(i,i+1);
                  end
                  if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     aa = aa + fit(i,i+1);
                  end
                  if (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;
               end
                  
                  
         if (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;


targ = fit;

%Start of main routine

n = size(prox,1);
acondiff = 1.0;
addcon = 0.0;
fit = prox(inperm,inperm);
proxave = zeros(n,n);
fmadd = zeros(n,n);

while (acondiff >= 1.0e-004)
   
   for i = 1:n
      for j = 1:n
         if(i ~= j)
            fit(i,j) = prox(inperm(i),inperm(j)) + addcon;
         else
            fit(i,j) = 0.0;
         end
      end
   end
   
   addconpv = addcon;
   
   work = zeros(n*(n-1)*(n-2),1);
   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;
         end
            
         if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     fit(i,i+1) = fit(i,i+1) - work(indexll+1);
                  end
                  if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     aac = aac + targ(i,i+1);
                  end
                  if (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);
                  end
                  if ((i >= jone) & (i < jtwo))
                     aa = aa + fit(i,i+1);
                  end
                  if (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;
               end
                  
                  
         if (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



addcon = -sum(sum(prox(inperm,inperm) - fit))/(n*(n-1));
acondiff = abs(addcon - addconpv);

end
   
  aveprox = sum(sum(prox))/(n*(n-1));

for i = 1:n
   for j = 1:n
      if ( i ~= j)
         proxave(i,j) = aveprox;
         fmadd(i,j) = fit(i,j) - addcon;
      else
         proxave(i,j) = 0;
         fmadd(i,j) = 0.0;
      end
   end
end

diff = sum(sum((prox(inperm,inperm) - (fmadd)).^2));

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

vaf = 1 - (diff/denom);
 
   
%toc;       
                        
                        
                     
                     
                     
                     
