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;