function pvalSimes=pval_Simes(pval,u)
%input
%pval:  a Vxn matrix of contrast p-values
%u: the partial conjunction parameter 
%(to find voxels where at least u active)
pvalSimes = NaN*ones(size(pval,1),1);
n=size(pval,2);   
for v=1:size(pval,1)
    p=squeeze(pval(v,:));
    n = sum(~isnan(p));
    simesFactor = (n-u+1)./(1:(n-u+1));
    if n>(u-1)
        p=sort(p(~isnan(p)));
        pvalSimes(v) = min(p(u:n).*simesFactor);
    end
end;