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

% QUADASSIGN 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 randomly permutated NSTARTS times with the
% cross-product index repeatedly computed with respect to PROXB.  The
% raw cross-product index between PROXA and PROXB is RAWINDEX; the
% NSTARTS 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 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,expectation,variance,upvalue,sortvalues] =
% quadassign(proxa,proxb,nstarts)

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

mean = ((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)));

indexvals = zeros(nstarts,1);

upval = 0;

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

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

