function [fit, diff, coord] = linfit_tied(prox,inperm,tiedcoord)

%  LINFIT_TIED does a confirmatory fitting of a given 
%  unidimensional order using Dykstra's
%  (Kaczmarz's) iterative projection least-squares method. This
%  includes the possible imposition of tied coordinates.
%
%  syntax: [fit, diff, coord] = linfit_tied(prox,inperm,tiedcoord)
%
%  INPERM is the given order;
%  FIT is an $n \times n$ matrix that is fitted to 
%  PROX(INPERM,INPERM) with least-squares value DIFF;
%  COORD gives the ordered coordinates whose absolute
%  differences could be used to reconstruct FIT; TIEDCOORD
%  is the tied pattern of coordinates imposed (in order)
%  along the continuum (using the integers from 1 up to n
%  to indicate the tied positions).

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

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
   
    for jone = 1:(n-1)
        for jtwo = (jone + 1):n
            p1 = fit(jone,jtwo);
            fit(jone,jtwo) = fit(jone,jtwo) - work(index + 1);
            
            if(tiedcoord(jone) == tiedcoord(jtwo))
                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

   
  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);

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

diff = (.5)*diff;