function ModSpec(ID, Method, Code)

%Modality specificity analysis

SaveDir = [pwd '/SampleResult/'];
DataDir = [pwd '/'];

load([SaveDir 'RidgeResult_' Method '_' num2str(Code) '_' ID '.mat'], 'Result');
ccs_w2w = Result.ccs1; %Word2Word
ccs_f2f = Result.ccs2; %Form2Form

load([SaveDir 'RidgeResult_CrossModal_' Method '_' num2str(Code) '_' ID '.mat'], 'cmResult');
ccs_w2f = cmResult.ccs_w2f; %Word2Form
ccs_f2w = cmResult.ccs_f2w; %Form2Word

%Negative value  to zero
ccs_w2w(find(ccs_w2w<0)) = 0;
ccs_f2f(find(ccs_f2f<0)) = 0;
ccs_w2f(find(ccs_w2f<0)) = 0;
ccs_f2w(find(ccs_f2w<0)) = 0;
  
%Modality specificity
MS_w = ccs_w2w - ccs_f2w;
MS_f = ccs_f2f - ccs_w2f;

%SaveResult
msResult.MS_w = MS_w;
msResult.MS_f = MS_f;

%%Significant testing
%Random seed
rng(1234,'twister')
%Permutation number
PermN = 10000;
[rMS_w, rMS_f] = Calc_rMS(Result, PermN);


%%Significant values of MS_w
% Transform to P values for each voxel
PX = [];
for ii = 1:length(MS_w)
    x = find(rMS_w>MS_w(ii));
    px = length(x)/PermN;    
    PX(ii) = px;
end

%FDR correction
Q = 0.05;
[PXsorted PXind] = sort(PX, 'ascend');
FDRthr = Q*[1:length(MS_w)]/length(MS_w);
Diff = PXsorted - FDRthr;
% Maximum index among non-significant voxels (in ascending order)        
if min(Diff) < 0
    if max(find(Diff<0)) == length(PXind)
        thrInd = PXind(end);
    else
        thrInd = PXind(max(find(Diff<0))+ 1);
    end
else
    thrInd = PXind(1);
end
thr = MS_w(thrInd);
disp(['Threshold MS(Word) = ' num2str(thr)])
% Modality Invariance, only significant voxels (determined by concatenated model pred acc)
NSvoxels = find(MS_w<=thr);
MS_w(NSvoxels) = 0;

msResult.MS_w_fdr = MS_w;
msResult.fdrthr_w = thr; 
msResult.NSvoxels_w = NSvoxels;

%Mapping from 1d Data to 3d .nii data 
Y = zeros(prod(Result.datasize),1);
for ii=1:length(Result.tvoxels)
    Y(Result.tvoxels(ii))= MS_w(ii);
end
vol = reshape(Y,Result.datasize);
vol_perm = permute(vol, [2,1,3]);
V = MRIread(expInfo.RefEPI);
V.vol = vol_perm;
MRIwrite(V,[SaveDir 'ModSpec_Word_'  Method '_' num2str(Code) '_' ID '_FDR.nii']);


%%Significant values of MS_f
% Transform to P values for each voxel
PX = [];
for ii = 1:length(MS_f)
    x = find(rMS_f>MS_f(ii));
    px = length(x)/PermN;    
    PX(ii) = px;
end

%FDR correction
Q = 0.05;
[PXsorted PXind] = sort(PX, 'ascend');
FDRthr = Q*[1:length(MS_f)]/length(MS_f);
Diff = PXsorted - FDRthr;
% Maximum index among non-significant voxels (in ascending order)        
if min(Diff) < 0
    if max(find(Diff<0)) == length(PXind)
        thrInd = PXind(end);
    else
        thrInd = PXind(max(find(Diff<0))+ 1);
    end
else
    thrInd = PXind(1);
end
thr = MS_f(thrInd);
disp(['Threshold MS(Word) = ' num2str(thr)])
% Modality Invariance, only significant voxels (determined by concatenated model pred acc)
NSvoxels = find(MS_f<=thr);
MS_f(NSvoxels) = 0;

msResult.MS_f_fdr = MS_f;
msResult.fdrthr_f = thr; 
msResult.NSvoxels_f = NSvoxels;

%Mapping from 1d Data to 3d .nii data
RefEPI = [DataDir 'target_' ID '.nii'];
Y = zeros(prod(Result.datasize),1);
for ii=1:length(Result.tvoxels)
    Y(Result.tvoxels(ii))= MS_f(ii);
end
vol = reshape(Y,Result.datasize);
vol_perm = permute(vol, [2,1,3]);
V = MRIread(RefEPI);
V.vol = vol_perm;
MRIwrite(V,[ 'ModSpec_Form_'  Method '_' num2str(Code) '_' ID '_FDR.nii']);

save([SaveDir 'ModSpec_' Method '_' num2str(Code) '_' ID '.mat'], 'msResult', '-v7.3');

end



function [rMS_w, rMS_f] = Calc_rMS(Result, PermN)

    A = normrnd(0,1,size(Result.resp1,1),PermN);
    B = normrnd(0,1,size(Result.resp1,1),PermN);
    rccs_w2w = mvn_corr_TN(A,B);

    A = normrnd(0,1,size(Result.resp2,1),PermN);
    B = normrnd(0,1,size(Result.resp2,1),PermN);
    rccs_f2f = mvn_corr_TN(A,B);

    A = normrnd(0,1,size(Result.resp2,1),PermN);
    B = normrnd(0,1,size(Result.resp2,1),PermN);
    rccs_w2f = mvn_corr_TN(A,B);
    
    A = normrnd(0,1,size(Result.resp1,1),PermN);
    B = normrnd(0,1,size(Result.resp1,1),PermN);
    rccs_f2w = mvn_corr_TN(A,B);

    %Negative value  to zero
    rccs_w2w(find(rccs_w2w<0)) = 0;
    rccs_f2f(find(rccs_f2f<0)) = 0;
    rccs_w2f(find(rccs_w2f<0)) = 0;
    rccs_f2w(find(rccs_f2w<0)) = 0;

    %Modality specificity
    rMS_w = rccs_w2w - rccs_f2w;
    rMS_f = rccs_f2f - rccs_w2f;        

end