function [rawindex,expectation,variance,upvalue,sortvalues] = quadassign_exact(proxa,proxb)

% QUADASSIGN_EXACT compares the two $n \times n$ square (but not necessarily
% symmetric) input matrices PROXA and PROXB (both with main diagonals of
% zero) using a (quadratic assignment) permutation strategy.  The rows
% and columns of PROXA are permutated in all possible
% $n!$ times (complete enumeration) with the
% cross-product index repeatedly computed with respect to PROXB.  The
% raw cross-product index between PROXA and PROXB is RAWINDEX; the
% $n!$ cross-product index values are given (sorted in increasing order)
% in the vector SORTVALUES; EXPECTATION and VARIANCE are the exact mean and 
% variance of the permutation distribution assuming that all $n!$ permutations
% are equally-likely; UPVALUE is the proportion of the $n!$ indices that
% are as large or larger than RAWINDEX (i.e., UPVALUE is the upper-tail p-value for
% the observed index RAWINDEX).  Given the type of computational storage required,
% $n$ should be restricted to being less than or equal to, say, nine.
%
% syntax:  [rawindex,expectation,variance,upvalue,sortvalues] =
% quadassign_exact(proxa,proxb)
n = size(proxa,1);

bone = ((sum(sum(proxa)))^2)*((sum(sum(proxb)))^2);

btwo = (sum(sum(proxa.^2)))*(sum(sum(proxb.^2)));

bthree = (sum(sum(proxa.*(proxa'))))*(sum(sum(proxb.*(proxb'))));

bfour = (sum(((sum(proxa)).^2) - sum(proxa.^2)))*(sum(((sum(proxb)).^2) - sum(proxb.^2)));

bfive = ((sum(proxa))*((sum(proxa'))') - sum(sum(proxa.*(proxa'))))* ...
    (sum(proxb)*((sum(proxb'))') - sum(sum(proxb.*(proxb'))));

bsix = (sum(((sum(proxa')).^2) - sum((proxa').^2)))* ...
    (sum(((sum(proxb')).^2) - sum((proxb').^2)));

bseven = (((sum(sum(proxa)))^2) - sum((sum(proxa')).^2) - 2*(sum(proxa)*((sum(proxa'))')) ...
    -sum((sum(proxa)).^2) + sum(sum(proxa.*(proxa'))) + sum(sum(proxa.^2)))* ...
     (((sum(sum(proxb)))^2) - sum((sum(proxb')).^2) - 2*(sum(proxb)*((sum(proxb'))')) ...
    -sum((sum(proxb)).^2) + sum(sum(proxb.*(proxb'))) + sum(sum(proxb.^2)));

rawindex = sum(sum(proxa.*proxb));

expectation = ((sum(sum(proxa)))*(sum(sum(proxb))))/(n*(n-1));

variance = (-bone/(n*n*(n-1)*(n-1))) + ((btwo + bthree)/(n*(n-1))) + ...
    ((bfour + 2*bfive + bsix)/(n*(n-1)*(n-2))) + (bseven/(n*(n-1)*(n-2)*(n-3)));

upval = 0;

permutations = perms(1:n);

nfactorial = size(permutations,1);
indexvals = zeros(nfactorial,1);

for i = 1:nfactorial
    
    perm = permutations(i,:);
    indexvals(i) = sum(sum((proxa(perm,perm)).*proxb));
    
    if (indexvals(i) >= rawindex)
        upval = upval + 1;
    end
end

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