%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Heller R., Golland Y., Malach R. and Benjamini Y., Conjunction Group   %
%Analysis: An alternative to mixed/random effect analysis,              %
%NeuroImage, Vol. 37(4), pp. 1178-1185, 2007.                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  


path(path,'C:\Users\ruheller\Conjunction\software-fMRI') %put here the directory where the matlab functions are

datadir='E:\fmri\data\malach_data\data070319\'  %put here the directory where your data is located in vmp format
filenames = {'DP_FP'; 'DP_HP'; 'DP_OP'}
cd 'E:\fmri\data\malach_data\temp'  %put here your working directory, where you want all the work to be saved.

%calculate a p-value map for all the t-maps in filenames
for i=1:length(filenames)
   tdata =readVMP(datadir,char(filenames(i)));
    tdata(tdata==0) = NaN;
    pmap =  1-normcdf(tdata);
    save(char(strcat('pmap_',filenames(i))),'pmap');
end    

pmap_names = {char(strcat('pmap_',filenames(1))); char(strcat('pmap_',filenames(2))); char(strcat('pmap_',filenames(3)))}%put here the maps to combine
u=[1; 2; 3] %put here the parameter u for screening (you can specify more than one such parameter)
pval_Combining(pmap_names,u,'Simes','pooledmap_DP');

%apply the FDR on each pooled p-value map
isActive = zeros(size(pooledmap,1), length(u));
unames = {'u1'; 'u2'; 'u3'}
for i=1:length(u)
    load(char(strcat('pooledmap_DP_',unames(i))))
    pID = fdr(pooledmap,0.05); 
    if ~isempty(pID)
        isActive(find(pooledmap<=pID),i) = 1;
    end
end
save('isActive_DP', 'isActive');

sum(isActive,1)%for u=1, 108301 voxels detected; for u=2, 12477 voxels detected; for u=3, 3317 voxels detected. 

%to control the OFDR if the screening hypothesis is u=1: 
load('pooledmap_DP_u1')
isActiveOFDR = zeros(size(pooledmap,1), 1);
pID = fdr(pooledmap,0.05); 
isActiveOFDR(find(pooledmap<=pID),1) = 1;
R = sum(isActive(:,1));
alpha = R*0.05/size(pooledmap,1)
for i=2:length(u)
    load(char(strcat('pooledmap_DP_',unames(i))))
    isActiveOFDR(find(pooledmap<=alpha),1) = i;  
end

tabulate(isActiveOFDR(isActiveOFDR>0))% for u=1, 

