function [rawindex,sampmean,sampvariance,upvalue,ncompare,sortvalues] = quartassign(proxa,proxb,nstarts)

% QUARTASSIGN compares the two $n \times n$ square (but not necessarily
% symmetric) input matrices PROXA and PROXB (both with main diagonals of
% zero) using a (quartic assignment) permutation strategy.  The rows
% and columns of PROXA are randomly permuted NSTARTS times with the
% four-place cross-product index repeatedly computed with respect to PROXB.  The
% raw four-place cross-product index (which is twice the number of consistent
% quadruples minus the inconsistent quadruples) between PROXA and PROXB is
% RAWINDEX; this is given explictly by 
% $\sum_{ i \ne j,k \ne l} sign(a_{kl} - a_{ij}) \times sign(b_{kl} - b_{ij}$
% NCOMPARE is the maximum RAWINDEX could be taking into account the tied
% proximity values in PROXA and PROXB; the NSTARTS cross-product index values
% are given (sorted in increasing order) in the vector SORTVALUES;
% SAMPMEAN and SAMPVARIANCE are the sample mean and variance of the NSTARTS
% values in SORTVALUES; UPVALUE is the proportion of the NSTARTS indices that
% are as large or larger than RAWINDEX (i.e., UPVALUE is the upper-tail p-value for
% the observed index RAWINDEX).
%
% syntax:  [rawindex,sampmean,sampvariance,upvalue,ncompare,sortvalues] =
% quartassign(proxa,proxb,nstarts)

n = size(proxa,1);

afourplace = zeros(n,n,n,n);
bfourplace = zeros(n,n,n,n);

for i = 1:n
    for j = 1:n
        for k = 1:n
            for l = 1:n
                if ( (i ~= j) & (k ~= l) )
                  afourplace(i,j,k,l) = sign(proxa(k,l) - proxa(i,j));
                  bfourplace(i,j,k,l) = sign(proxb(k,l) - proxb(i,j));
                end
            end
        end
     end
 end


rawindex = sum(sum(sum(sum(afourplace.*bfourplace))));

ncompare = sum(sum(sum(sum(abs(afourplace.*bfourplace)))));

indexvals = zeros(nstarts,1);

upval = 0;

for i = 1:nstarts
    
    perm = randperm(n);
    indexvals(i) = sum(sum(sum(sum((afourplace(perm,perm,perm,perm).*bfourplace)))));
    
    if (indexvals(i) >= rawindex)
        upval = upval + 1;
    end
end

sortvalues = sort(indexvals)';
upvalue = upval/nstarts;

sampmean = mean(sortvalues);
sampvariance = (std(sortvalues,1))^2;