function [outpermone,outpermtwo,coordone,coordtwo,fitone,fittwo,addconone,addcontwo,dev] = ...
biscallp(prox,inpermone,inpermtwo,nopt)


%BISCALLP carries out a bidimensional scaling of a symmetric proximity
%  matrix using iterative linear programming.
%  PROX is the input proximity matrix (with a zero main diagonal and a 
%  dissimilarity interpretation);
%  INPERMONE is the input beginning permutation for the first dimension
%  (a permuation of the first $n$ integers); INPERMTWO is the input beginning
%  permutation for the second dimension;
%  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 LP
%  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;
%  DEV 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 = zeros(n,n);
fittwo = zeros(n,n);


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


proxone = prox;
proxtwo = zeros(n,n);
proxmed = zeros(n,n);

nch2 = n*(n-1)/2;
proxvec = zeros(nch2,1);

index = 0;

for i = 2:n
    for j = 1:(i-1)
        index = index + 1;
        proxvec(index) = prox(i,j);
    end
end
medprox = median(proxvec);



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


devdiff = 1.0;
dev = 0.0;

while (devdiff >= 1.0e-005)
   
devprev = dev;

if (nopt == 1)
    [outperm,coordone,devone,fitone,addconone] = ...
        uniscallpac(proxone,outpermone);
    outpermone = outperm;
end

if (nopt == 0)
    [fitone,devone,coordone,addconone,exitflag] = linfitl1ac(proxone,outpermone);
end

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



if (nopt == 1)
    [outperm,coordtwo,devtwo,fittwo,addcontwo] = ...
        uniscallpac(proxtwo,outpermtwo);
    outpermtwo = outperm;
end

if (nopt == 0)
    [fittwo,devtwo,coordtwo,addcontwo,exitflag] = linfitl1ac(proxtwo,outpermtwo);
end

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);
            
        else
            proxthree(outpermtwo(i),outpermtwo(j)) = 0.0;
        end
    end
end

              
diff = sum(sum(abs(proxthree)));
denom = sum(sum(abs(prox-proxmed)));
dev = 1 - (diff/denom)

devdiff = abs(dev-devprev);

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

toc