function [fit, vaf, coord, addcon, eqfit, eqvaf, eqaddcon] = ...
           eqspace_linfitac(prox,inperm)

%  EQSPACE_LINFITAC does a confirmatory fitting of a given unidimensional order
%  using the Dykstra--Kaczmarz iterative projection 
%  least-squares method, but differing from linfit.m in 
%  including the estimation of an additive constant.  Also
%  an equally-spaced confirmatory fitting alternative is carried out.
%
%  syntax: [fit, vaf, coord, addcon, eqfit, eqvaf, eqaddcon] = ...
%             eqspace_linfitac(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 VAF;
%  COORD gives the ordered coordinates whose absolute differences
%  could be used to reconstruct FIT; ADDCON is the estimated
%  additive constant that can be interpreted as being added to PROX.
%  The equally-spaced output alternatives are prefixed with an EQ.


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

while (acondiff >= 1.0e-005)
   
   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;

while (cr >= 1.0e-006)
   
   cr = 0.0;
   index = 0;
   
   for jone = 1:(n-2)
      for jtwo = (jone+1):(n-1)
         for jthree = (jtwo+1):n
            
            p1 = fit(jone,jtwo);
            p2 = fit(jone,jthree);
            p3 = fit(jtwo,jthree);
            
            fit(jone,jtwo) = fit(jone,jtwo) - work(index+1);
            fit(jone,jthree) = fit(jone,jthree) - work(index+2);
            fit(jtwo,jthree) = fit(jtwo,jthree) - work(index+3);
            
            del = (fit(jone,jthree) - fit(jone,jtwo) - ...
               fit(jtwo,jthree))/3.0;
            
            fit(jone,jthree) = fit(jone,jthree) - del;
            fit(jone,jtwo) = fit(jone,jtwo) + del;
            fit(jtwo,jthree) = fit(jtwo,jthree) + del;
            
            work(index+1) = del;
            work(index+2) = -del;
            work(index+3) = del;
            
            index = index + 3;
            
            cr = cr + abs(p1-fit(jone,jtwo)) + abs(p2-fit(jone,jthree)) ...
               + abs(p3 - fit(jtwo,jthree));
            
         end
      end
   end
   
   for jone = 1:(n-1)
      for jtwo = (jone+1):n
         
         p1 = fit(jone,jtwo);
         
         fit(jone,jtwo) = fit(jone,jtwo) - work(index+1);
         
         if(fit(jone,jtwo) < 0.0)
            
            work(index+1) = -fit(jone,jtwo);
            fit(jone,jtwo) = 0.0;
            
        else
            
            work(index+1) = 0.0;
            
         end
         
         index = index + 1;
         
         cr= cr + abs(p1-fit(jone,jtwo));
      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);

coord = zeros(n,1);

for i = 1:n
    for j = 1:n
        if(i>j)
            coord(i) = coord(i) + fit(i,j);
        end
        if(i<j)
            coord(i) = coord(i) - fit(i,j);
        end
    end
    
    coord(i) = coord(i)/n;
end

exog = zeros(n,n);

for i = 1:n
    for j = 1:n
        exog(i,j) = abs(i-j);
    end
end

xsum = 0.0;
xsumsq = 0.0;
ysum = 0.0;
ysumsq = 0.0;
xysum = 0.0;

for i = 1:(n-1)
    for j = (i+1):n
        
        xsum = xsum + exog(i,j);
        xsumsq = xsumsq + exog(i,j)*exog(i,j);
        ysum = ysum + prox(inperm(i),inperm(j));
        ysumsq = ysumsq + prox(inperm(i),inperm(j))*prox(inperm(i),inperm(j));
        xysum = xysum + exog(i,j)*prox(inperm(i),inperm(j));
        
    end
end

nct = n*(n-1)/2;

byx = (xysum - ((1/nct)*xsum*ysum))/(xsumsq - ((1/nct)*xsum*xsum));

eqaddcon = ((1/nct)*ysum) - (byx*(1/nct)*xsum);

eqfit = zeros(n,n);

for i = 1:n
    for j = 1:n
        
        if(i ~= j)
            
            eqfit(i,j) = eqaddcon + byx*exog(i,j);
            
                     
        else
            eqfit(i,j) = 0;
           
           
        end
    end
end

eqfit;

eqdiff = sum(sum((prox(inperm,inperm) - eqfit).^2));

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

eqvaf = 1 - (eqdiff/eqdenom);

for i = 1:n
    for j = 1:n
        
        if(i ~= j)
            
            eqfit(i,j) = byx*exog(i,j);
            
        else
            eqfit(i,j) = 0;
        end
    end
end

eqaddcon = -eqaddcon;
