function p=PCDVp(Gr1,Gr2,nEval)
% p=PCDVp(Gr1,Gr2,nEval)    or     p=PCDVp(Data,n1,nEval)
% In the first form, Gr1 and Gr2 are data from two independent groups of cases and have identical dimensions, except possibly for the number of
% levels of the last dimension, which correspond to the number of cases in each group
% In the second form, all the data are in Data (with cases as last index) and n1 indicates that the first n1 cases from Data belong to one group, 
% the remainder belonging to the other group
% if Data or Gr1 and Gr2 have more than two dimensions, they are made two dimensional by merging all but their last
% dimension
% A projection score on the mean difference between groups is first calculated for the original group membership and
% then on nEval-1 random permutations of group memberships (keeping each group size constant)
% p is the probability of the data under the null hypothesis that group membership is irrelevant
% It is obtained as the rank of the projection score on the original data divided by nEval
% Any score equality with the score on the original data is given rang priority over the original one (avoiding inflated
% type 1 error, especially when the number of cases is small so that fewer than nEval different sign permutations exist)
% As with any permutation test, different runs on the same data may produce different probability estimates
% reference:
% Achim, A.  Statistical detection of between groups differences in event-related potentials.  Clinical Neurophysiology, 2001, 112: 1023-1034.

a=size(Gr1);
n1=a(end);
k=length(a);
if length(Gr2)>1,
    b=size(Gr2);
    if length(b)~=k | any(a(1:k-1)~=b(1:k-1)), error('The two matrices passed to PCDVp are not compatible.'); end 
    if k>2,  %if 3 dimensions, transform to 2 dimensions
        D=[reshape(Gr1,prod(a(1:k-1)),n1) reshape(Gr2,prod(b(1:k-1)),b(k))];
    else
        D=[Gr1 Gr2];    
    end
else
    if length(a)>2, D=reshape(Gr1,prod(a(1:k-1)),n1); else D=Gr1; end
    n1=Gr2;
end

nn=size(D,2);
s0=project(D,n1);   % score on actual group membership
p=1;
for n=2:nEval,
   if project(D(:,randperm(nn)),n1)>=s0, p=p+1; end
end
p=p/nEval;



function crit=project(X,n);
% crit is related (motonic positive) to a Student t for independent groups, which is OK since only ranks count
% (computations accelerated by squaring the numerator while keeping its sign, rather than taking sqrt in denominator
% and it is useless to devide all scores by same constant)
nn=size(X,2);
n1=n+1;
m1=mean(X(:,1:n),2);
m2=mean(X(:,n1:nn),2);
d=m1-m2;
s=X'*d;
m1=mean(s(1:n));
m2=mean(s(n1:nn));
s(1:n)=s(1:n)-m1;
s(n1:nn)=s(n1:nn)-m2;
crit=(m1-m2).^2/(s'*s); 
if m2>m1, crit=-crit; end