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

% CUBEASSIGN compares the two $n \times n$ square (but not necessarily
% symmetric) input matrices PROXA and PROXB (both with main diagonals of
% zero) using a (cubic assignment) permutation strategy.  The rows
% and columns of PROXA are randomly permuted NSTARTS times with the
% triad cross-product index repeatedly computed with respect to PROXB.  The
% raw triad cross-product index (which is twice the number of consistent
% row triads minus the inconsistent row triads) between PROXA and PROXB is
% RAWINDEX; this is given explictly by 
% $\sum_{distinct i,j,k} sign(a_{ik} - a_{ij}) \times sign(b_{ik} - 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).

n = size(proxa,1);

athreeplace = zeros(n,n,n);
bthreeplace = zeros(n,n,n);

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

rawindex = sum(sum(sum(athreeplace.*bthreeplace)));

ncompare = sum(sum(sum(abs(athreeplace.*bthreeplace))));

indexvals = zeros(nstarts,1);

upval = 0;

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

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

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