function [fit,vaf,rowperm,colperm,addcon,coord] = linfittmac(proxtm,inperm)

%  LINFITTMAC does a confirmatory two-mode fitting of a given 
%  unidimensional ordering of the row and column objects of
%  a two-mode proximity matrix PROXTM using Dykstra's (Kaczmarz's)
%  iterative projection least-squares method;
%  it differs from linfittm.m  by including the estimation of an
%  additive constant.
%
%  syntax: [fit,vaf,rowperm,colperm,addcon,coord] = ...
%    linfittmac(proxtm,inperm)
%
%  INPERM is the given ordering of the row and column objects
%  together; FIT is an nrow (number of rows) by ncol (number
%  of columns) matrix  of absolute coordinate differences that
%  is fitted to PROXTM(ROWPERM,COLPERM) with VAF being the
%  variance-accounted-for.  ROWPERM and COLPERM are the row and
%  column object orderings derived from INPERM.  ADDCON is the
%  estimated additive constant that can be interpreted as being
%  added to PROXTM (or alternatively subtracted
%  from the fitted matrix FIT).  The nrow + ncol coordinates
%  (ordered with the smallest
%  set at a value of zero) are given in COORD.


[nrow ncol] = size(proxtm);
n = nrow + ncol;
tempfit = zeros(n,n);
fit = zeros(nrow,ncol);
mfit = zeros(n,n);

for irow = 1:nrow
   for icol = 1:ncol
      
      tempfit(irow,icol+nrow) = proxtm(irow,icol);
      tempfit(icol+nrow,irow) = proxtm(irow,icol);
      
   end
end

acondiff = 1.0;
addcon = 0.0;

while (acondiff >= 1.0e-005)

  for i =1:n
      for j=1:n
         
         if(i ~= j)
            mfit(i,j) = tempfit(inperm(i),inperm(j)) + addcon;
         else
            mfit(i,j) = 0.0;
         end
      end
   end

cr = 1.0;
work = zeros(n*(n-1)*(n-2)*(n-3),1);
addconpv = addcon;

while(cr >= 1.0e-005)
   
   cr = 0.0;
   indexll=0;
   
   for ione = 1:(n-3)
      for itwo = (ione+1):(n-2)
         for ithree = (itwo+1):(n-1)
            for ifour = (ithree+1):n
               
               if((((inperm(ione) <= nrow) && (inperm(itwo) <= nrow)) ...
                     && ((inperm(ithree) > nrow) && (inperm(ifour) ...
                     > nrow))) || ...
                     (((inperm(ione) > nrow) && (inperm(itwo) > nrow)) ...
                     && ((inperm(ithree) <= nrow) && (inperm(ifour) ...
                     <= nrow))))
                  
                  p1 = mfit(ione,ithree);
                  p2 = mfit(ione,ifour);
                  p3 = mfit(itwo,ithree);
                  p4 = mfit(itwo,ifour);
                  
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+1);
                  mfit(ione,ifour) = mfit(ione,ifour) - work(indexll+2);
                  mfit(itwo,ithree) = mfit(itwo,ithree) - work(indexll+3);
                  mfit(itwo,ifour) = mfit(itwo,ifour) - work(indexll+4);
                  
                  del = (mfit(itwo,ithree) + mfit(ione,ifour) - ...
                     mfit(ione,ithree) - mfit(itwo,ifour))/4.0;
                  
                  mfit(itwo,ithree) = mfit(itwo,ithree) - del;
                  mfit(ione,ifour) = mfit(ione,ifour) - del;
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(itwo,ifour) = mfit(itwo,ifour) + del;
                  
                  work(indexll+1) = del;
                  work(indexll+2) = -del;
                  work(indexll+3) = -del;
                  work(indexll+4) = del;
                  
                  cr = cr + abs(p1 - mfit(ione,ithree)) + ...
                     abs(p2 - mfit(ione,ifour)) + ...
                     abs(p3 - mfit(itwo,ithree)) + ...
                     abs(p4 - mfit(itwo,ifour));
                  
                  indexll = indexll + 4;
                  
                  
                elseif((((inperm(ione) <= nrow) && (inperm(itwo) > nrow)) ...
                  && ((inperm(ithree) <= nrow) && (inperm(ifour) ...
                  > nrow))) || ...
                  (((inperm(ione) > nrow) && (inperm(itwo) <= nrow)) ...
                  && ((inperm(ithree) > nrow) && (inperm(ifour) ...
                  <= nrow))))
                  
                  p1 = mfit(ione,itwo);
                  p2 = mfit(ione,ifour);
                  p3 = mfit(itwo,ithree);
                  p4 = mfit(ithree,ifour);
                  
                  mfit(ione,itwo) = mfit(ione,itwo) - work(indexll+1);
                  mfit(ione,ifour) = mfit(ione,ifour) - work(indexll+2);
                  mfit(itwo,ithree) = mfit(itwo,ithree) - work(indexll+3);
                  mfit(ithree,ifour) = mfit(ithree,ifour) - work(indexll+4);
                  
                  del = (mfit(ione,itwo) + mfit(itwo,ithree) + ...
                     mfit(ithree,ifour) - mfit(ione,ifour))/4.0;
                  
                  mfit(itwo,ithree) = mfit(itwo,ithree) - del;
                  mfit(ione,ifour) = mfit(ione,ifour) + del;
                  mfit(ione,itwo) = mfit(ione,itwo) - del;
                  mfit(ithree,ifour) = mfit(ithree,ifour) - del;
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = del;
                  work(indexll+3) = -del;
                  work(indexll+4) = -del;
                  
                  cr = cr + abs(p1 - mfit(ione,itwo)) + ...
                     abs(p2 - mfit(ione,ifour)) + ...
                     abs(p3 - mfit(itwo,ithree)) + ...
                     abs(p4 - mfit(ithree,ifour));
                  
                  indexll = indexll + 4;
                  
                elseif((((inperm(ione) <= nrow) && (inperm(itwo) > nrow)) ...
                  && ((inperm(ithree) > nrow) && (inperm(ifour) ...
                  <= nrow))) || ...
                  (((inperm(ione) > nrow) && (inperm(itwo) <= nrow)) ...
                  && ((inperm(ithree) <= nrow) && (inperm(ifour) ...
                  > nrow))))
                  
                  p1 = mfit(ione,itwo);
                  p2 = mfit(ione,ithree);
                  p3 = mfit(itwo,ifour);
                  p4 = mfit(ithree,ifour);
                  
                  mfit(ione,itwo) = mfit(ione,itwo) - work(indexll+1);
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+2);
                  mfit(itwo,ifour) = mfit(itwo,ifour) - work(indexll+3);
                  mfit(ithree,ifour) = mfit(ithree,ifour) - work(indexll+4);
                  
                  del = (mfit(ione,itwo) + mfit(itwo,ifour) - ...
                     mfit(ione,ithree) - mfit(ithree,ifour))/4.0;
                  
                  mfit(ione,itwo) = mfit(ione,itwo) - del;
                  mfit(itwo,ifour) = mfit(itwo,ifour) - del;
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(ithree,ifour) = mfit(ithree,ifour) + del;
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = del;
                  work(indexll+3) = -del;
                  work(indexll+4) = del;
                  
                  cr = cr + abs(p1 - mfit(ione,itwo)) + ...
                     abs(p2 - mfit(ione,ithree)) + ...
                     abs(p3 - mfit(itwo,ifour)) + ...
                     abs(p4 - mfit(ithree,ifour));
                  
                  indexll = indexll + 4;
                  
               end
               
            end
         end
      end
   end

   for ione = 1:(n-2)
      for itwo = (ione+1):(n-1)
         for ithree = (itwo+1):n
                          
               if((((inperm(ione) <= nrow) && (inperm(itwo) <= nrow)) ...
                     && ((inperm(ithree) > nrow) )) || ...
                     (((inperm(ione) > nrow) && (inperm(itwo) > nrow)) ...
                     && ((inperm(ithree) <= nrow) )))
                  
                  p1 = mfit(ione,ithree);
                  p2 = mfit(itwo,ithree);
                                    
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+1);
                  mfit(itwo,ithree) = mfit(itwo,ithree) - work(indexll+2);
                 
                  if (mfit(ione,ithree) >= mfit(itwo,ithree))
                      del = 0;
                  else
                      del = (mfit(itwo,ithree) - mfit(ione,ithree))/2;
                  end
                  
                               
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(itwo,ithree) = mfit(itwo,ithree) - del;
                  
                  
                  work(indexll+1) = del;
                  work(indexll+2) = -del;
                  
                  
                  cr = cr + abs(p1 - mfit(ione,ithree)) + ...
                     abs(p2 - mfit(itwo,ithree));
                     
                 indexll = indexll + 2;
                  
                
                elseif((((inperm(ione) <= nrow) && (inperm(itwo) > nrow)) ...
                  && ((inperm(ithree) > nrow) )) || ...
                  (((inperm(ione) > nrow) && (inperm(itwo) <= nrow)) ...
                  && ((inperm(ithree) <= nrow) )))
                  
                  p1 = mfit(ione,itwo);
                  p2 = mfit(ione,ithree);
                  
                  
                  mfit(ione,itwo) = mfit(ione,itwo) - work(indexll+1);
                  mfit(ione,ithree) = mfit(ione,ithree) - work(indexll+2);
                  
                  if (mfit(ione,ithree) >= mfit(ione,itwo))
                      del = 0;
                  else
                      del = (mfit(ione,itwo) - mfit(ione,ithree))/2;
                  end
                  
                               
                  mfit(ione,ithree) = mfit(ione,ithree) + del;
                  mfit(ione,itwo) = mfit(ione,itwo) - del;
                  
                  
                  work(indexll+1) = -del;
                  work(indexll+2) = del;
                  
                  
                  cr = cr + abs(p1 - mfit(ione,itwo)) + ...
                     abs(p2 - mfit(ione,ithree));
                     
                 indexll = indexll + 2;
                  
                  
                  
               end
               
            end
         end
      end
      
   

   
   
   for ione = 1:(n-1)
      for itwo = (ione+1):n
         
         if(((inperm(ione) <= nrow) && (inperm(itwo) > nrow)) || ...
               ((inperm(ione) > nrow) && (inperm(itwo) <= nrow)))
            
            p1 = mfit(ione,itwo);
            
            mfit(ione,itwo) = mfit(ione,itwo) - work( indexll+1);
            
            if(mfit(ione,itwo) < 0.0)
               
               work(indexll+1) = -mfit(ione,itwo);
               mfit(ione,itwo) = 0.0;
            end
            
            indexll = indexll + 1;
            
            cr = cr + abs(p1 - mfit(ione,itwo));
         end
      end
      
   end
end

   
   
   for i = 1:(n-1)
      for j = (i+1):n
         
         mfit(j,i) = mfit(i,j);
         
      end
   end
   
   rowidx = 0;
   
   for i = 1:n
      
      if(inperm(i) <= nrow)
         
         rowidx = rowidx + 1;
         
         colidx = 0;
         
         for j = 1:n
            
            if(inperm(j) > nrow)
               
               colidx = colidx + 1;
               fit(rowidx,colidx) = mfit(i,j);
            end
         end
      end
   end
   
   rowperm = zeros(nrow,1);
   colperm = zeros(ncol,1);
   rowidx = 0;
   colidx = 0;
   
   for i = 1:n
      
      if(inperm(i) <= nrow)
         
         rowidx = rowidx + 1;
         rowperm(rowidx) = inperm(i);
         
      elseif(inperm(i) > nrow)
         
         colidx = colidx + 1;
         colperm(colidx) = inperm(i) - nrow;
         
      end
   end
   
    
   addcon = -sum(sum(proxtm(rowperm,colperm) - fit))/(nrow*ncol);

   acondiff = abs(addcon - addconpv);

   aveprox = sum(sum(proxtm))/(nrow*ncol);
   
   diff = sum(sum((proxtm(rowperm,colperm) - (fit-addcon)).^2));
   denom = sum(sum((proxtm(rowperm,colperm) - aveprox).^2));
   
   vaf = 1 - (diff/denom);
  
end  
      
    coord = zeros(n,1);
   
   if (inperm(1) <= nrow)
       
       for i = 2:n
           if(inperm(i) > nrow)
               coord(i) = mfit(1,i);
           end
       end
       
       for i = 2:(n-1)
           for j = (i+1):n
               if( (inperm(i) <= nrow) && (inperm(j) > nrow ) )
                   coord(i) = coord(j) - mfit(i,j);
               end
               if( (inperm(i) > nrow) && (inperm(j) <= nrow) )
                   coord(j) = coord(i) + mfit(i,j);
               end
           end
       end
   end
      
   
  if (inperm(1) >  nrow)
       
       for i = 2:n
           if(inperm(i) <= nrow)
               coord(i) = mfit(1,i);
           end
       end
       
       for i = 2:(n-1)
           for j = (i+1):n
               if( (inperm(i) <= nrow) && (inperm(j) > nrow ) )
                   coord(j) = coord(i) + mfit(i,j);
               end
               if( (inperm(i) > nrow) && (inperm(j) <= nrow) )
                   coord(i) = coord(j) - mfit(i,j);
               end
           end
       end
   end
      
   
       
  
   
       