function [z,zn,p]=jMantelTest(A,B,nPerms)
% 
% Conduct normalised Mantel hypothesis test for the null association of
% two square symmetric matrices
%
% IN
%   A, B: square symmetric matrices of real numbers
%   nPerms: scalar, number of permutations
%
% OUT
%   z: Mantel test statistic (cross-product of individual matrix entries)
%   zn: normalised Mantel test statistic (Smouse et al)
%   p: p-value for the null hypothesis of no association between the
%   matrices, based on the normalised statistic
%
% REFERENCES 
% - N. Mantel, The Detection of Disease Clustering and a Generalized
%   Regression Approach, Cancer Research (27), pp. 209-220, 1967
% - P. E. Smouse et al., Multiple Regression and Correlation Extensions
%   of the Mantel Test of Matrix Correspondence, Systematic Zoology(V35N4), 
%   pp.627-632, 1986
% - H. E. Daniels, The relation between measures of correlation in the
%   universe of sample permutations, Biometrika, pp.129-135, 1944
%
% NOTE
% - The null hypothesis is not well defined in the Mantel paper, the
%   "null hypothesis of no clustering" is what comes closest.
%
% VERSION
% v1.0 oct 2013 Jonas Richiardi
% http://www.stanford.edu/~richiard/software.html
% - Initial release
% v1.1 oct 2013 Jonas Richiardi & Will Shirer
% - Handle negative relationships: matrices can be similar even though
%   they anticovary
% v1.1.1 29 oct 2013 JR
% - fixed bug in p-value computation

% minimal sanity check
if any(isnan(A(:))) || any(isnan(B(:)))
    warning('jMantelTest:nans',['data contains nans, replacing with 0s' ...
        ', please check your data and mistrust results']);
    A(isnan(A))=0; B(isnan(B))=0;
end
% this test fails if match is not exact.
if ~(isequal(triu(A),tril(A)') && isequal(triu(B),tril(B)'))
    warning('jMantelTest:sym',['Both matrices are not symmetric (up ' ...
        'to machine precision), please check your data and mistrust ' ...
        'results']);
end

Av=jUpperTriMatToVec(A,1);     % vector of upper tri mat of A
Bv=jUpperTriMatToVec(B,1);     % vector of upper tri mat of B
zpn=zeros(nPerms,1);           % permutation values of the normalised test statistic

% compute test stat values
z=Av'*Bv;                     % Mantel test stat
zn=corr(Av,Bv);               % Normalised Mantel test stats

if zn<0
    warning('jMantelTest:neg','Matrices have a negative relationship.');
end

for i=1:nPerms
    % permute rows & cols of matrix B
    rndIdx=randperm(size(B,1));
    Bvp=jUpperTriMatToVec(B(rndIdx,rndIdx),1);
    zpn(i)=corr(Av,Bvp);
end
p=(1+sum(abs(zpn)>=abs(zn)))/(nPerms+1);