#!/usr/bin/env python
# coding: utf-8

import json
import csv
import PredysUtil8ans as pu
import numpy as np
import os
from sklearn.svm import SVC
from sklearn.svm import LinearSVC
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import accuracy_score
import nibabel as nib
from nibabel import load
import math    
import csv
import glob
import random
import pickle
import copy
import scipy
from scipy import signal
import scipy.stats
from scipy import ndimage
from statsmodels.stats.multitest import fdrcorrection
    


def IntraModalDec(TaskName, MethodName='Test', TargTime=[0,1,2,3,4], bdelay=3, Rad=5, RunFlag=2, svmType=1):

#Input:
# MethodName ... output label
# TargTime ... Target time points used
# bdelay ... Delay of bold signal from onset
# Rad ... Radius of searchlight sphere
# RunFlag ... 1: two runs are used (if available) / 2: Only the first run is used
# svmType ... 1: SVC with RBF kernel / 2: LinearSVC 

    if (TaskName == 'Words'):
        SubList = [1,3,4,5,6,7,9,10,11,12,13,14,17,18,19,20,25,28,29,30,31,32,37,38,40,41,43,44,45,46,48,50,55,58,59,61,64,67,68,70,72] 
    else:
        SubList = [1,2,3,4,5,6,7,9,10,11,12,13,14,15,17,18,19,20,24,25,27,28,29,30,31,32,37,38,40,41,43,44,45,46,48,50,54,55,58,59,61,64,67,68,70,72] 

    CurrDir = os.getcwd()
    RoiDir = CurrDir + '/WFU_AALAtlas'
    SaveDir = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName    
    PreprocDataDir = CurrDir + '/DetrendData_8ans'        
    RoiName = 'GrayMatter'
    MinVoxN = 10 #Minimum voxel numbers for a sphere
    VolumeSize = [79, 95, 39]
    TotalVoxelSize = np.prod(VolumeSize)
    if RunFlag == 1:
        UseRun = [1, 2]
    else:
        UseRun = [1]

    #Load Whole brain data
    RoiImg = load(os.path.join(RoiDir,RoiName + '.nii'))
    RoiVol = RoiImg.get_fdata()
    Roi_x, Roi_y, Roi_z = np.nonzero(RoiVol)
    RoiLen = len(Roi_x)

    #Indices of voxels included in Sphere ROIs    
    VoxelIndex = [np.unravel_index(ii,VolumeSize) for ii in range(TotalVoxelSize)]
    TargVoxel = [ii for ii in range(TotalVoxelSize) if RoiVol[VoxelIndex[ii]] > 0]
    VoxelIndex_targ = [VoxelIndex[TargVoxel[tt]] for tt in range(len(TargVoxel))]

    #Load preprocessed brain response data
    tRespVal01 = []
    tRespVal02 = []
    for SubNum in SubList:
        SubName = 'sub-{:02d}'.format(SubNum)
        print('Loading ' + SubName +' ...')
        fname = 'RespData_bd{:d}_{}_NoAdapt_{}_Run{:02d}.npy'
        if (TaskName == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName, SubName, RunNum), TargTime) for RunNum in [2]]
        else:
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName, SubName, RunNum), TargTime) for RunNum in UseRun]
        T01_List = [T01_List[tt] for tt in range(len(T01_List)) if not len(T01_List[tt])==0]
        T01 = np.concatenate(T01_List, axis=0)
        tRespVal01.append(T01)

        fname = 'RespData_bd{:d}_{}_Adapt_{}_Run{:02d}.npy' 
        if (TaskName == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName, SubName, RunNum), TargTime) for RunNum in [2]]
        else:        
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName, SubName, RunNum), TargTime) for RunNum in UseRun]
        T02_List = [T02_List[tt] for tt in range(len(T02_List)) if not len(T02_List[tt])==0]
        T02 = np.concatenate(T02_List, axis=0)
        tRespVal02.append(T02)

    #Load searchlight sphere data
    with open('SearchSphereSet_Rad{}.pickle'.format(Rad), mode='rb') as f:
        SphereSet = pickle.load(f) 

    #Select svm type
    if svmType == 1:
        model = make_pipeline(StandardScaler(), SVC(gamma='auto')) 
    elif svmType == 2:
        model = make_pipeline(StandardScaler(), LinearSVC(random_state=0))

    #Prepare saving directory
    if not (os.path.exists(SaveDir)):
        os.mkdir(SaveDir)

    #Searchlight analysis
    tData = np.zeros([RoiLen,len(SubList)])
    for rr in range(RoiLen): #Loop for voxels in whole brain
        print('Calculating searchligth at voxel #' + str(rr) )
        #Get indices of target voxels included in Sphere ROI
        SpherePos = SphereSet[rr]
        vind = [VoxelIndex_targ.index(tuple(SpherePos[vv])) for vv in range(len(SpherePos))]
        if len(vind) < MinVoxN: #Skip if number of voxels in the shpere is very small
            tAcc = np.NaN
        else:
            tAcc = np.zeros(len(SubList))
            #Leave-one-subject-out cross-validation
            for TestSub in range(len(SubList)):
                TrainingSub = list(range(len(SubList)))
                del TrainingSub[TestSub]

                #Concatenate training subjects data
                Val01_List = [tRespVal01[ss][:,vind] for ss in TrainingSub]
                X01 = np.concatenate(Val01_List, axis=0)
                Val02_List = [tRespVal02[ss][:,vind] for ss in TrainingSub]
                X02 = np.concatenate(Val02_List, axis=0)
                X = np.concatenate([X01,X02],axis=0)   #Concatenate data from two conditions
                X = np.nan_to_num(X) #Replace NaN to 0
                Y = np.concatenate([np.zeros(X01.shape[0]), np.ones(X02.shape[0])])

                #Concatenate test subjects data
                X01_test = tRespVal01[TestSub][:,vind]
                X02_test = tRespVal02[TestSub][:,vind]
                X_test = np.concatenate([X01_test,X02_test],axis=0)
                X_test = np.nan_to_num(X_test)
                Y_test = np.concatenate([np.zeros(X01_test.shape[0]), np.ones(X02_test.shape[0])])

                #SVM model training
                model.fit(X, Y)
                #Model testing
                Y_pred = model.predict(X_test)
                tAcc[TestSub]  = accuracy_score(Y_pred,Y_test)

        print('Mean Accuracy: {}'.format(round(np.mean(tAcc),3)))  
        tData[rr,:] = tAcc
        
    fname = 'SearchLight_LOOCV_{}_{}'.format(TaskName, MethodName)
    np.save(os.path.join(SaveDir,fname),tData)
    



def CrossModalDec(TaskName1, TaskName2, MethodName='Test', TargTime=[0,1,2,3,4], bdelay=3, Rad=5, RunFlag=2, svmType=1):

#Input:
# MethodName ... output label
# TargTime ... target time points used
# bdelay ... delay of bold signal from onset
# Rad ... radius of searchlight sphere
# RunFlag ... 1: two runs are used (if available) / 2: Only the first run is used
# svmType ... 1: SVC with RBF kernel / 2: LinearSVC 

    SaveTaskName = TaskName1 + '2' + TaskName2
    if (TaskName1 == 'Words') or (TaskName2 == 'Words'):
        SubList = [1,3,4,5,6,7,9,10,11,12,13,14,17,18,19,20,25,28,29,30,31,32,37,38,40,41,43,44,45,46,48,50,55,58,59,61,64,67,68,70,72] 
    else:
        SubList = [1,2,3,4,5,6,7,9,10,11,12,13,14,15,17,18,19,20,24,25,27,28,29,30,31,32,37,38,40,41,43,44,45,46,48,50,54,55,58,59,61,64,67,68,70,72] 

    CurrDir = os.getcwd()
    RoiDir = CurrDir + '/WFU_AALAtlas'
    SaveDir = CurrDir + '/Result_8ans/Result_SearchLight_' + SaveTaskName    
    PreprocDataDir = CurrDir + '/DetrendData_8ans'        
    RoiName = 'GrayMatter'
    MinVoxN = 10 #Minimum voxel numbers for a sphere
    VolumeSize = [79, 95, 39]
    TotalVoxelSize = np.prod(VolumeSize)
    if RunFlag == 1:
        UseRun = [1, 2]
    else:
        UseRun = [1]

    #Load Whole brain data
    RoiImg = load(os.path.join(RoiDir,RoiName + '.nii'))
    RoiVol = RoiImg.get_fdata()
    Roi_x, Roi_y, Roi_z = np.nonzero(RoiVol)
    RoiLen = len(Roi_x)

    #Indices of voxels included in Sphere ROIs    
    VoxelIndex = [np.unravel_index(ii,VolumeSize) for ii in range(TotalVoxelSize)]
    TargVoxel = [ii for ii in range(TotalVoxelSize) if RoiVol[VoxelIndex[ii]] > 0]
    VoxelIndex_targ = [VoxelIndex[TargVoxel[tt]] for tt in range(len(TargVoxel))]      

    #Load training data
    tRespVal01_trn = []
    tRespVal02_trn = []
    for SubNum in SubList:
        SubName = 'sub-{:02d}'.format(SubNum)
        print('Loading ' + SubName +' ...')
        fname = 'RespData_bd{:d}_{}_NoAdapt_{}_Run{:02d}.npy'         
        if (TaskName1 == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName1, SubName, RunNum), TargTime) for RunNum in [2]]
        else:
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName1, SubName, RunNum), TargTime) for RunNum in UseRun]
        T01_List = [T01_List[tt] for tt in range(len(T01_List)) if not len(T01_List[tt])==0]
        T01 = np.concatenate(T01_List, axis=0)
        tRespVal01_trn.append(T01)

        fname = 'RespData_bd{:d}_{}_Adapt_{}_Run{:02d}.npy' 
        if (TaskName1 == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName1, SubName, RunNum), TargTime) for RunNum in [2]]
        else: 
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName1, SubName, RunNum), TargTime) for RunNum in UseRun]
        T02_List = [T02_List[tt] for tt in range(len(T02_List)) if not len(T02_List[tt])==0]
        T02 = np.concatenate(T02_List, axis=0)
        tRespVal02_trn.append(T02)

    #Load test data
    tRespVal01_test = []
    tRespVal02_test = []
    for SubNum in SubList:
        SubName = 'sub-{:02d}'.format(SubNum)
        print('Loading ' + SubName +' ...')
        fname = 'RespData_bd{:d}_{}_NoAdapt_{}_Run{:02d}.npy'         
        if (TaskName2 == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName2, SubName, RunNum), TargTime) for RunNum in [2]]
        else:
            T01_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName2, SubName, RunNum), TargTime) for RunNum in UseRun]
        T01_List = [T01_List[tt] for tt in range(len(T01_List)) if not len(T01_List[tt])==0]
        T01 = np.concatenate(T01_List, axis=0)
        tRespVal01_test.append(T01)

        fname = 'RespData_bd{:d}_{}_Adapt_{}_Run{:02d}.npy' 
        if (TaskName2 == 'Words') and (RunFlag == 2) and (SubNum == 70): #For this subject, only 2nd run exists
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName2, SubName, RunNum), TargTime) for RunNum in [2]]
        else: 
            T02_List = [pu.CalcTrialRespVal(PreprocDataDir, fname.format(bdelay, TaskName2, SubName, RunNum), TargTime) for RunNum in UseRun]
        T02_List = [T02_List[tt] for tt in range(len(T02_List)) if not len(T02_List[tt])==0]
        T02 = np.concatenate(T02_List, axis=0)
        tRespVal02_test.append(T02)

    #Load searchlight sphere data
    with open('SearchSphereSet_Rad{}.pickle'.format(Rad), mode='rb') as f:
        SphereSet = pickle.load(f) 

    #Select svm type
    if svmType == 1:
        model = make_pipeline(StandardScaler(), SVC(gamma='auto')) 
    elif svmType == 2:
        model = make_pipeline(StandardScaler(), LinearSVC(random_state=0))

    #Prepare saving directory
    if not (os.path.exists(SaveDir)):
        os.mkdir(SaveDir)
        
    #Searchlight analysis
    tData = np.zeros([RoiLen,len(SubList)])
    for rr in range(RoiLen): #Loop for voxels in whole brain
        print('Calculating searchligth at voxel #' + str(rr) )
        #Get indices of target voxels included in Sphere ROI
        SpherePos = SphereSet[rr]
        vind = [VoxelIndex_targ.index(tuple(SpherePos[vv])) for vv in range(len(SpherePos))]
        if len(vind) < MinVoxN: #Skip if number of voxels in the shpere is very small
            tAcc = np.NaN
        else:
            tAcc = np.zeros(len(SubList))
            #Leave-one-subject-out cross-validation
            for TestSub in range(len(SubList)):
                TrainingSub = list(range(len(SubList)))
                del TrainingSub[TestSub]

                #Concatenate training subjects data
                Val01_List = [tRespVal01_trn[ss][:,vind] for ss in TrainingSub]
                X01 = np.concatenate(Val01_List, axis=0)
                Val02_List = [tRespVal02_trn[ss][:,vind] for ss in TrainingSub]
                X02 = np.concatenate(Val02_List, axis=0)
                X = np.concatenate([X01,X02],axis=0)   #Concatenate data from two conditions
                X = np.nan_to_num(X) #Replace NaN to 0
                Y = np.concatenate([np.zeros(X01.shape[0]), np.ones(X02.shape[0])])

                #Concatenate test subjects data
                X01_test = tRespVal01_test[TestSub][:,vind]
                X02_test = tRespVal02_test[TestSub][:,vind]
                X_test = np.concatenate([X01_test,X02_test],axis=0)
                X_test = np.nan_to_num(X_test)
                Y_test = np.concatenate([np.zeros(X01_test.shape[0]), np.ones(X02_test.shape[0])])

                #SVM model training
                model.fit(X, Y)
                #Model testing
                Y_pred = model.predict(X_test)
                tAcc[TestSub]  = accuracy_score(Y_pred,Y_test)

        print('Mean Accuracy: {}'.format(round(np.mean(tAcc),3)))  
        tData[rr,:] = tAcc
        
    fname = 'SearchLight_LOOCV_{}_{}'.format(SaveTaskName, MethodName)
    np.save(os.path.join(SaveDir, fname), tData)



def SaveNifti_PermTest(TaskName, MethodName='Test', PeakThr=0.005, ClusterThr=0.05, N_Rand=50000):

    #N_Rand = 50000 #Number of permutation
    #PeakThr = 0.005 #Peak level threshold
    #ClusterThr = 0.05 #Cluster level threshold

    CurrDir = os.getcwd()
    SaveDir = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName + '/'
    fname = os.path.join(SaveDir, 'SearchLight_LOOCV_' + TaskName + '_' + MethodName + '.npy')
    SaveName = 'SearchLight_LOOCV_{}_{}_PermTest_LogPval.nii'.format(TaskName, MethodName)
    SaveName_FDR = 'SearchLight_LOOCV_{}_{}_PermTest_Peak{}_FDR.nii'.format(TaskName, MethodName, str(PeakThr)[2:])
    
    #Load decoding result
    tData = np.load(fname)
    
    #Sign permutation test
    pData, mData_rn = pu.SignPerm(tData, N_Rand) 

    #Transform 1d p value result into 3d MNI space
    RoiName = 'GrayMatter.nii'
    RoiDir = CurrDir + '/WFU_AALAtlas'
    RoiImg = nib.load(os.path.join(RoiDir, RoiName))
    RoiVol = RoiImg.get_fdata()

    #Null distribution of cluster size
    random.seed(1234)     
    N_Rand_Cluster = round(N_Rand / 100)
    tClusterSize = [pu.CalcClusterSize(tData, mData_rn, RoiVol, PeakThr, rr) for rr in range(N_Rand_Cluster)]
    ClusterSizeDist = np.concatenate(tClusterSize)

    #Calculate actual cluster size
    pData_vol = np.ones(np.shape(RoiVol))
    pData_vol[np.nonzero(RoiVol)] = pData
    cluster_map, n_clusters = ndimage.label(pData_vol < PeakThr)
    TargClusterSize = np.bincount(cluster_map.ravel())[1:]

    #Transform cluster size data to p values
    ClusterPval = np.zeros(np.shape(TargClusterSize)[0])
    for cc in range(np.shape(TargClusterSize)[0]):
        ClusterPval[cc] = np.count_nonzero(ClusterSizeDist >= TargClusterSize[cc]) / np.shape(ClusterSizeDist)[0]

    #FDR correction
    fdr_mask, p_fdr = fdrcorrection(ClusterPval, alpha=ClusterThr)

    #Masking p value with fdr clustere correction
    Mask = np.zeros(np.shape(pData_vol))
    for mm in np.where(fdr_mask)[0]:
        Cluster = np.where(cluster_map==mm+1)
        for vv in range(np.shape(Cluster)[1]):
            Mask[Cluster[0][vv],Cluster[1][vv],Cluster[2][vv]] = 1

    #Transform p values to -log(p)
    pData_vol[pData_vol==0] = 1 / N_Rand # Correction in case of pvalue = 0
    pData_vol = abs(np.log10(pData_vol))
    pData_vol_fdr = pData_vol*Mask

    #Save data
    SaveImg = nib.Nifti1Image(pData_vol, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName)
    SaveImg = nib.Nifti1Image(pData_vol_fdr, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName_FDR)
    
    return

    
def SaveNifti_PermTest_Conj(TaskName1, TaskName2, MethodName='Test', PeakThr=0.005, ClusterThr=0.05, N_Rand=50000):

    #Conjunction analysis with sign permutation

    #N_Rand = 50000 #Number of permutation
    #PeakThr = 0.005 #Peak level threshold
    #ClusterThr = 0.05 #Cluster level threshold

    CurrDir = os.getcwd()
    LoadDir1 = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName1 + '/'
    fname1 = os.path.join(LoadDir1, 'SearchLight_LOOCV_' + TaskName1 + '_' + MethodName + '.npy')
    LoadDir2 = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName2 + '/'
    fname2 = os.path.join(LoadDir2, 'SearchLight_LOOCV_' + TaskName2 + '_' + MethodName + '.npy')   
    SaveDir = CurrDir + '/Result_8ans/Result_SearchLight_Conjunction/'
    SaveName = 'SearchLight_LOOCV_Conjunction_{}_{}_{}_PermTest_LogPval.nii'.format(TaskName1, TaskName2, MethodName)
    SaveName_FDR = 'SearchLight_LOOCV_Conjunction_{}_{}_{}_PermTest_Peak{}_FDR.nii'.format(TaskName1, TaskName2, MethodName, str(PeakThr)[2:])

    #Prepare saving directri
    if not (os.path.exists(SaveDir)):
        os.mkdir(SaveDir)

    #Load decoding result
    tData1 = np.load(fname1)
    tData2 = np.load(fname2)

    #Sign permutation test
    pData1, mData_rn1 = pu.SignPerm(tData1, N_Rand) 
    pData2, mData_rn2 = pu.SignPerm(tData2, N_Rand)

    #Get max p value of 2 conditions
    pData_max = np.max(np.stack([pData1, pData2]),0)

    #Transform 1d p value result into 3d MNI space
    RoiName = 'GrayMatter.nii'
    RoiDir = CurrDir + '/WFU_AALAtlas'
    RoiImg = nib.load(os.path.join(RoiDir, RoiName))
    RoiVol = RoiImg.get_fdata()
    
    #Null distribution of cluster size
    random.seed(1234)     
    N_Rand_Cluster = round(N_Rand / 100)
    tClusterSize = [pu.CalcClusterSize_Conj(tData1, tData2, mData_rn1, mData_rn2, RoiVol, PeakThr, rr) for rr in range(N_Rand_Cluster)]
    ClusterSizeDist = np.concatenate(tClusterSize)

    #Calculate actual cluster size
    pData_vol = np.ones(np.shape(RoiVol))
    pData_vol[np.nonzero(RoiVol)] = pData_max
    cluster_map, n_clusters = ndimage.label(pData_vol < PeakThr)
    TargClusterSize = np.bincount(cluster_map.ravel())[1:]

    #Transform cluster size data to p values
    ClusterPval = np.zeros(np.shape(TargClusterSize)[0])
    for cc in range(np.shape(TargClusterSize)[0]):
        ClusterPval[cc] = np.count_nonzero(ClusterSizeDist >= TargClusterSize[cc]) / np.shape(ClusterSizeDist)[0]

    #FDR correction
    fdr_mask, p_fdr = fdrcorrection(ClusterPval, alpha=ClusterThr)

    #Masking p value with fdr cluster correction
    Mask = np.zeros(np.shape(pData_vol))
    for mm in np.where(fdr_mask)[0]:
        Cluster = np.where(cluster_map == mm+1)
        for vv in range(np.shape(Cluster)[1]):
            Mask[Cluster[0][vv], Cluster[1][vv], Cluster[2][vv]] = 1

    #Transform p values to -log(p)
    pData_vol[pData_vol==0] = 1 / N_Rand # Correction in case of pvalue = 0
    pData_vol = abs(np.log10(pData_vol))
    pData_vol_fdr = pData_vol * Mask

    #Save data
    SaveImg = nib.Nifti1Image(pData_vol, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName)
    SaveImg = nib.Nifti1Image(pData_vol_fdr, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName_FDR)
    
    return



def SaveNifti_Paired_PermTest(TaskName, MethodName='Test', PeakThr=0.005, ClusterThr=0.05, N_Rand=50000):

    #Two-sample sign permutation test (8yo vs 5yo)
    
    #For 8 yo data, only one of two runs is used
    #N_Rand = 50000 #Number of permutation
    #PeakThr = 0.005 #Peak level threshold
    #ClusterThr = 0.05 #Cluster level threshold

    CurrDir = os.getcwd()
    SaveDir1 = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName + '/'
    SaveDir2 = CurrDir + '/Result_5ans/Result_SearchLight_' + TaskName + '/'    
    fname1 = os.path.join(SaveDir1, 'SearchLight_LOOCV_' + TaskName + '_' + MethodName + '.npy')
    fname2 = os.path.join(SaveDir2, 'SearchLight_LOOCV_' + TaskName + '_' + MethodName + '.npy')    
    SaveName = 'SearchLight_LOOCV_Paired_8yo-5yo_{}_{}_PermTest_LogPval.nii'.format(TaskName, MethodName)
    SaveName_FDR = 'SearchLight_LOOCV_Paired_8yo-5yo_{}_{}_PermTest_Peak{}_FDR.nii'.format(TaskName, MethodName, str(PeakThr)[2:])
    
    #Load decoding result
    tData1 = np.load(fname1)
    tData2 = np.load(fname2)

    #Sign permutation test
    pData, mData_rn = pu.SignPerm_Paired(tData1, tData2, N_Rand) 

    #Transform 1d p value result into 3d MNI space
    RoiName = 'GrayMatter.nii'
    RoiDir = CurrDir + '/WFU_AALAtlas'
    RoiImg = nib.load(os.path.join(RoiDir, RoiName))
    RoiVol = RoiImg.get_fdata()
    
    #Null distribution of cluster size
    random.seed(1234) 
    N_Rand_Cluster = round(N_Rand / 100)
    tClusterSize = [pu.CalcClusterSize_Paired(tData1, tData2, mData_rn, RoiVol, PeakThr, rr) for rr in range(N_Rand_Cluster)]
    ClusterSizeDist = np.concatenate(tClusterSize)

    #Calculate actual cluster size
    pData_vol = np.ones(np.shape(RoiVol))
    pData_vol[np.nonzero(RoiVol)] = pData
    cluster_map, n_clusters = ndimage.label(pData_vol < PeakThr)
    TargClusterSize = np.bincount(cluster_map.ravel())[1:]

    #Transform cluster size data to p values
    ClusterPval = np.zeros(np.shape(TargClusterSize)[0])
    for cc in range(np.shape(TargClusterSize)[0]):
        ClusterPval[cc] = np.count_nonzero(ClusterSizeDist >= TargClusterSize[cc]) / np.shape(ClusterSizeDist)[0]

    #FDR correction
    fdr_mask, p_fdr = fdrcorrection(ClusterPval, alpha=ClusterThr)

    #Masking p value with fdr clustere correction
    Mask = np.zeros(np.shape(pData_vol))
    for mm in np.where(fdr_mask)[0]:
        Cluster = np.where(cluster_map == mm+1)
        for vv in range(np.shape(Cluster)[1]):
            Mask[Cluster[0][vv], Cluster[1][vv], Cluster[2][vv]] = 1

    #Transform p values to -log(p)
    pData_vol[pData_vol==0] = 1 / N_Rand # Correction in case of pvalue = 0
    pData_vol = abs(np.log10(pData_vol))
    pData_vol_fdr = pData_vol * Mask

    #Save data
    SaveImg = nib.Nifti1Image(pData_vol, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir1 + SaveName)
    SaveImg = nib.Nifti1Image(pData_vol_fdr, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir1 + SaveName_FDR)
    
    return


def SaveNifti_Paired_PermTest_Conj(TaskName1, TaskName2, MethodName='Test', PeakThr=0.005, ClusterThr=0.05, N_Rand=50000):

    #Two-sample sign permutation test (8yo vs 5yo) with Conjunction analysis

    #N_Rand = 50000 #Number of permutation
    #PeakThr = 0.005 #Peak level threshold
    #ClusterThr = 0.05 #Cluster level threshold

    CurrDir = os.getcwd()
    SaveDir = CurrDir + '/Result_8ans/Result_SearchLight_Conjunction/'
    LoadDir1A = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName1 + '/'
    LoadDir1B = CurrDir + '/Result_5ans/Result_SearchLight_' + TaskName1 + '/'
    LoadDir2A = CurrDir + '/Result_8ans/Result_SearchLight_' + TaskName2 + '/'
    LoadDir2B = CurrDir + '/Result_5ans/Result_SearchLight_' + TaskName2 + '/'      
    fname1A = os.path.join(LoadDir1A, 'SearchLight_LOOCV_' + TaskName1 + '_' + MethodName + '.npy')
    fname1B = os.path.join(LoadDir1B, 'SearchLight_LOOCV_' + TaskName1 + '_' + MethodName + '.npy')
    fname2A = os.path.join(LoadDir2A, 'SearchLight_LOOCV_' + TaskName2 + '_' + MethodName + '.npy')
    fname2B = os.path.join(LoadDir2B, 'SearchLight_LOOCV_' + TaskName2 + '_' + MethodName + '.npy')       
    SaveName = 'SearchLight_LOOCV_Paired_8yo-5yo_Conjunction_{}_{}_{}_PermTest_LogPval.nii'.format(TaskName1, TaskName2, MethodName)
    SaveName_FDR = 'SearchLight_LOOCV_Paired_8yo-5yo_Conjunction_{}_{}_{}_PermTest_Peak{}_FDR.nii'.format(TaskName1, TaskName2, MethodName, str(PeakThr)[2:])
    
    #Load decoding result
    tData1A = np.load(fname1A)
    tData1B = np.load(fname1B)
    tData2A = np.load(fname2A)
    tData2B = np.load(fname2B)

    #Sign permutation test
    pData1, mData_rn1 = pu.SignPerm_Paired(tData1A, tData1B, N_Rand)
    pData2, mData_rn2 = pu.SignPerm_Paired(tData2A, tData2B, N_Rand) 

    #Get max p value of 2 conditions
    pData_max = np.max(np.stack([pData1, pData2]), 0)

    #Transform 1d p value result into 3d MNI space
    RoiName = 'GrayMatter.nii'
    RoiDir = CurrDir + '/WFU_AALAtlas'
    RoiImg = nib.load(os.path.join(RoiDir, RoiName))
    RoiVol = RoiImg.get_fdata()
    
    #Null distribution of cluster size
    random.seed(1234) 
    N_Rand_Cluster = round(N_Rand / 100)
    tClusterSize = [pu.CalcClusterSize_Paired_Conj(tData1A, tData1B, tData2A, tData2B, mData_rn1, mData_rn2, RoiVol, PeakThr, rr) for rr in range(N_Rand_Cluster)]
    ClusterSizeDist = np.concatenate(tClusterSize)

    #Calculate actual cluster size
    pData_vol = np.ones(np.shape(RoiVol))
    pData_vol[np.nonzero(RoiVol)] = pData_max
    cluster_map, n_clusters = ndimage.label(pData_vol < PeakThr)
    TargClusterSize = np.bincount(cluster_map.ravel())[1:]

    #Transform cluster size data to p values
    ClusterPval = np.zeros(np.shape(TargClusterSize)[0])
    for cc in range(np.shape(TargClusterSize)[0]):
        ClusterPval[cc] = np.count_nonzero(ClusterSizeDist >= TargClusterSize[cc]) / np.shape(ClusterSizeDist)[0]

    #FDR correction
    fdr_mask, p_fdr = fdrcorrection(ClusterPval, alpha=ClusterThr)

    #Masking p value with fdr clustere correction
    Mask = np.zeros(np.shape(pData_vol))
    for mm in np.where(fdr_mask)[0]:
        Cluster = np.where(cluster_map == mm+1)
        for vv in range(np.shape(Cluster)[1]):
            Mask[Cluster[0][vv], Cluster[1][vv], Cluster[2][vv]] = 1

    #Transform p values to -log(p)
    pData_vol[pData_vol==0] = 1 / N_Rand # Correction in case of pvalue = 0
    pData_vol = abs(np.log10(pData_vol))
    pData_vol_fdr = pData_vol * Mask

    #Save data
    SaveImg = nib.Nifti1Image(pData_vol, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName)
    SaveImg = nib.Nifti1Image(pData_vol_fdr, RoiImg.affine, RoiImg.header)
    nib.save(SaveImg, SaveDir + SaveName_FDR)
    
    return