function [outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo, ...
    addconone,addcontwo,vaf] = ...
ms_biscalqa(prox,targone,targtwo,inpermone,inpermtwo,kblock,nopt)

%  BISCALQA carries out a bidimensional scaling of a symmetric 
%  proximity matrix using iterative quadratic assignment.
%
%  syntax: [outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo,...
%    addconone,addcontwo,vaf] = ...
%     biscalqa(prox,targone,targtwo,inpermone,inpermtwo,kblock,nopt)
%
%  PROX is the input proximity matrix (with a zero main diagonal 
%  and a dissimilarity interpretation);
%  TARGONE is the input target matrix for the first dimension 
%  (usually with a zero main diagonal and a dissimilarity 
%  interpretation representing equally spaced locations along 
%  a continuum); TARGTWO is the input target
%  matrix for the second dimension;
%  INPERMONE is the input beginning permutation for the first 
%  dimension (a permutation of the first $n$ integers);
%  INPERMTWO is the input beginning
%  permutation for the second dimension;
%  the insertion and rotation routines use from 1 to KBLOCK
%  (which is less than or equal to $n-1$) consecutive objects in 
%  the permutation defining the row and column orders of the data
%  matrix. NOPT controls the confirmatory or exploratory fitting 
%  of the unidimensional scales; a value of NOPT = 0 will fit in a
%  confirmatory manner the two scales
%  indicated by INPERMONE and INPERMTWO;
%  a value of NOPT = 1 uses iterative QA
%  to locate the better permutations to fit;
%  OUTPERMONE is the final object permutation for the
%  first dimension; OUTPERMTWO is the final object permutation
%  for the second dimension;
%  COORDONE is the set of first dimension coordinates
%  in ascending order; COORDTWO is the set of second dimension
%  coordinates in ascending order;
%  ADDCONONE is the additive constant for the first
%  dimensional model; ADDCONTWO is the additive constant for 
%  the second dimensional model;
%  VAF is the variance-accounted-for in PROX by 
%  the bidimensional scaling.

%tic;
n = size(prox,1);
outpermone = inpermone;
outpermtwo = inpermtwo;
coordone = zeros(n,1);
coordtwo = zeros(n,1);
fitone = targone;
fittwo = targtwo;

addconone = 0.0;
addcontwo = 0.0;
fitonedim = zeros(n,n);
fittwodim = zeros(n,n);


proxone = prox;
proxtwo = zeros(n,n);
proxaveorig = zeros(n,n);
proxave = zeros(n,n);
proxthree = zeros(n,n);
aveproxorig = sum(sum(prox))/(n*(n-1));

for i = 1:n
    for j = 1:n
        if (i ~= j)
            proxaveorig(i,j) = aveproxorig;
        else
            proxaveorig(i,j) = 0.0;
            
        end
    end
end


vafdiff = 1.0;
vaf = 0.0;

while (vafdiff >= 1.0e-005)
   
vafprev = vaf;

if (nopt == 1)

begindexone = sum(sum(proxone(outpermone,outpermone).*fitone));

    nchange = 1;
   
while (nchange == 1)
   
   nchange=0;
   
   for k = 1:(n-1)
      for j = (k+1):n
                   
         intrperm = outpermone;
             
         intrperm(k) = outpermone(j);
         intrperm(j) = outpermone(k);
         
         tryindex = sum(sum(proxone(intrperm,intrperm).*fitone));
            
         if(tryindex > (begindexone + 1.0e-008))
            nchange = 1;
            begindexone = tryindex;
            outpermone = intrperm;
        end
         
      end
   end
   
   
   for k = 1:kblock
      for insertpt = 1:(n+1)
         for nlimlow = 1:(n+1-k)
            
            intrperm = outpermone;
            
            if (nlimlow > insertpt)
               
               jtwo = 0;
               for j = insertpt:(insertpt+k-1)
                  intrperm(j) =outpermone(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = (insertpt+k):(nlimlow+k-1);
                  intrperm(j) = outpermone(insertpt+jone);
                  jone = jone + 1;
               end
            end
               
            if ((nlimlow+k) < insertpt)
               
               jtwo = 0;
               for j = (insertpt-k):(insertpt-1)
                  intrperm(j) = outpermone(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = nlimlow:(insertpt-k-1)
                  intrperm(j) = outpermone(nlimlow+k+jone);
                  jone = jone + 1;
               end
               
            
               
            end
            
            tryindex = sum(sum(proxone(intrperm,intrperm).*fitone));
            
            if (tryindex > (begindexone + 1.0e-008))
               nchange = 1;
               begindexone = tryindex;
               outpermone = intrperm;
           end
            
         end
      end
   end
   

if (kblock > 1)
    
   for k = 2:kblock
      for nlimlow = 1:(n+1-k)
                   
         intrperm = outpermone;
            
         for j = 1:k
            intrperm(nlimlow+j-1) = outpermone(nlimlow+k-j);
         end
         
         tryindex = sum(sum(proxone(intrperm,intrperm).*fitone));
            
         if (tryindex > (begindexone + 1.0e-008))
            nchange = 1;
            begindexone = tryindex;
            outpermone = intrperm;
        end
         
      end
   end
end
end
end


%[fitone,vafone,coordone,addconone] = linfitac(proxone,outpermone);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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) = proxone(outpermone(i),outpermone(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(proxone(outpermone,outpermone) - fit))/(n*(n-1));

acondiff = abs(addcon - addconpv);

end

   
  aveprox = sum(sum(proxone))/(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((proxone(outpermone,outpermone) - (fmadd)).^2));

denom = sum(sum((proxone(outpermone,outpermone) - 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

fitone = fit;
coordone = coord;
vafone = vaf;
addconone = addcon;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fitonedim = fitone;

for i = 1:n
    for j = 1:n
        
        if(i ~= j)
            
            proxtwo(outpermone(i),outpermone(j)) = proxone(outpermone(i),outpermone(j)) - ...
                fitonedim(i,j) + addconone;
            
        else
            proxtwo(outpermone(i),outpermone(j)) = 0.0;
        end
    end
end
for i = 1:n
    for j = 1:n
        if (i ~= j)
            proxtwo(outpermtwo(i),outpermtwo(j)) = ...
                proxtwo(outpermtwo(i),outpermtwo(j)) + fittwodim(i,j) - addcontwo;
        else
            
        end
    end
end

if (nopt == 1)
    
  begindextwo = sum(sum(proxtwo(outpermtwo,outpermtwo).*fittwo));

    nchange = 1;
   
while (nchange == 1)
   
   nchange=0;
   
   for k = 1:(n-1)
      for j = (k+1):n
                   
         intrperm = outpermtwo;
             
         intrperm(k) = outpermtwo(j);
         intrperm(j) = outpermtwo(k);
         
         tryindex = sum(sum(proxtwo(intrperm,intrperm).*fittwo));
            
         if (tryindex > (begindextwo + 1.0e-008))
            nchange = 1;
            begindextwo = tryindex;
            outpermtwo = intrperm;
        end
         
      end
   end
   
   
   for k = 1:kblock
      for insertpt = 1:(n+1)
         for nlimlow = 1:(n+1-k)
            
            intrperm = outpermtwo;
            
            if (nlimlow > insertpt)
               
               jtwo = 0;
               for j = insertpt:(insertpt+k-1)
                  intrperm(j) =outpermtwo(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = (insertpt+k):(nlimlow+k-1);
                  intrperm(j) = outpermtwo(insertpt+jone);
                  jone = jone + 1;
               end
            end
               
            if ((nlimlow+k) < insertpt)
               
               jtwo = 0;
               for j = (insertpt-k):(insertpt-1)
                  intrperm(j) = outpermtwo(nlimlow+jtwo);
                  jtwo = jtwo + 1;
               end
               
               jone = 0;
               for j = nlimlow:(insertpt-k-1)
                  intrperm(j) = outpermtwo(nlimlow+k+jone);
                  jone = jone + 1;
               end
               
            
               
            end
            
            tryindex = sum(sum(proxtwo(intrperm,intrperm).*fittwo));
            
            if (tryindex > (begindextwo + 1.0e-008))
               nchange = 1;
               begindextwo = tryindex;
               outpermtwo = intrperm;
           end
            
         end
      end
   end
   

if (kblock > 1)
    
   for k = 2:kblock
      for nlimlow = 1:(n+1-k)
                   
         intrperm = outpermtwo;
            
         for j = 1:k
            intrperm(nlimlow+j-1) = outpermtwo(nlimlow+k-j);
         end
         
         tryindex = sum(sum(proxtwo(intrperm,intrperm).*fittwo));
            
         if (tryindex > (begindextwo + 1.0e-008))
            nchange = 1;
            begindextwo = tryindex;
            outpermtwo = intrperm;
        end
         
      end
   end
end
end
end

%[fittwo,vaftwo,coordtwo,addcontwo] = linfitac(proxtwo,outpermtwo);

%%%%%%%%%%%%%%%%%%%%%
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) = proxtwo(outpermtwo(i),outpermtwo(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(proxtwo(outpermtwo,outpermtwo) - fit))/(n*(n-1));

acondiff = abs(addcon - addconpv);

end

   
  aveprox = sum(sum(proxtwo))/(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((proxtwo(outpermtwo,outpermtwo) - (fmadd)).^2));

denom = sum(sum((proxtwo(outpermtwo,outpermtwo) - 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
%%%%%%%%%%%%%%%%%%%%%%

fittwo = fit;
coordtwo = coord;
vaftwo = vaf;
addcontwo = addcon;


fittwodim = fittwo;

for i = 1:n
    for j = 1:n
        
        if (i ~= j)
            
            proxthree(outpermtwo(i),outpermtwo(j)) = proxtwo(outpermtwo(i),outpermtwo(j)) - ...
                fittwodim(i,j) + addcontwo;
            
        else
            proxthree(outpermtwo(i),outpermtwo(j)) = 0.0;
        end
    end
end

              
diff = sum(sum(proxthree.^2));
denom = sum(sum((prox-proxaveorig).^2));
vaf = 1 - (diff/denom);

vafdiff = abs(vaf-vafprev);

for i = 1:n
    for j = 1:n
        
        if (i ~= j)
                    
        proxone(outpermone(i),outpermone(j)) = fitone(i,j) - addconone + ...
            proxthree(outpermone(i),outpermone(j));
        end
    
        if (i == j)
            
            proxone(outpermone(i),outpermone(j)) = 0.0;
            
        end
    end
end
            
end

%toc