Nipype-FSL fMRI Analysis

#load fsl module
import module
await module.load('fsl/6.0.4') #Original pipeline: FSL 6.0.4, nipype 1.6.1.
await module.list()

import os
os.environ["FSLOUTPUTTYPE"]="NIFTI_GZ"

# nipype 1.11.0's MultipleRegressDesign formats (N,1) numpy arrays with %e/%d,
# which numpy>=2 rejects ("only 0-dimensional arrays can be converted to Python
# scalars"). There is no numpy<2 wheel for Python 3.13, so rather than downgrade
# numpy (a slow source build), keep numpy 2.x and fix the two affected lines.
import sys, subprocess, glob, pathlib
subprocess.run([sys.executable, "-m", "pip", "install",
                "pandas==2.3.3", "nibabel==5.3.3", "nilearn==0.13.1", "nipype==1.11.0"],
               check=False)
for _m in map(pathlib.Path, glob.glob(f"{sys.prefix}/lib/python*/site-packages/nipype/interfaces/fsl/model.py")):
    _m.write_text(_m.read_text().replace("val for val in convals", "val for val in convals.ravel()"))
    print("Patched nipype MultipleRegressDesign for numpy>=2:", _m)

from __future__ import print_function
from __future__ import division
from builtins import str
from builtins import range

from glob import glob

from nipype import Function
from nipype.interfaces import fsl, utility as util, io as nio
import nipype.pipeline.engine as pe  # pypeline engine
import nipype.algorithms.modelgen as model  # model generation
import nipype.algorithms.rapidart as ra  # artifact detection
import matplotlib.pyplot as plt

import nilearn
from nilearn.image import math_img
from nilearn import surface, plotting, datasets
import nibabel as nib
import numpy as np

# check output type and fsl version
print(fsl.Info.output_type())
print(fsl.Info.version())

# Input: Set the BIDS path
data_dir_first = os.path.join(os.getcwd(), 'ds000102')
# Output: Set path where nipype will store stepwise results
exp_dir_first = os.path.join(os.getcwd(), 'output_level1')

# Input: Set path of first level outputs
data_dir_second = os.path.join(os.getcwd(), 'output_level1/level1_results/')
# Output: Set path where nipype will store stepwise results
exp_dir_second = os.path.join(os.getcwd(), 'output_level2/')

# Input: Set path of second level outputs
data_dir_third = os.path.join(os.getcwd(), 'output_level2/level2_results/')
# Output: Set path where nipype will store stepwise results
exp_dir_third = os.path.join(os.getcwd(), 'output_level3/')

PATTERN = "sub-0*"

# Use Nilearn's MNI152 template to keep the notebook self-contained and portable
mni = datasets.load_mni152_template(resolution=2)

# Corresponding brain mask
mask = datasets.load_mni152_brain_mask(resolution=2)

# Brain-extracted MNI template
mni_img = math_img("img1 * img2", img1=mni, img2=mask)

# Write to disk (portable, CI-safe)
mni_brain = os.path.abspath("MNI152_T1_2mm_brain_nilearn.nii.gz")
nib.save(mni_img, mni_brain)

plotting.plot_anat(mni_brain)

!datalad install https://github.com/OpenNeuroDatasets/ds000102.git
!cd ds000102 && datalad get $PATTERN

wf_first = pe.Workflow(name='level1', base_dir=exp_dir_first)
wf_first.config["execution"]["crashfile_format"] = "txt"

# get subject_id list - we only do 9 subjects here so the notebook runs in a reasonable time
subj_list = [x.split('-')[-1] for x in glob(data_dir_first+"/"+PATTERN)]
subj_list.sort()
print(subj_list)

infosource_first = pe.Node(util.IdentityInterface(fields=["subject_id"]),
                  name="infosource")
infosource_first.iterables = [("subject_id", subj_list)]

dg_first = pe.Node(
    interface=nio.DataGrabber(
        infields=["subject_id","run_id"], outfields=["struct", "func", "events"]
    ),
    name="dg_first"
)

# Specify run_ids and return a sorted filelist to ensure we match files to correct runs/tasks
dg_first.inputs.run_id = [1,2]
dg_first.inputs.sort_filelist = True
dg_first.inputs.template = "*"
dg_first.inputs.base_directory = data_dir_first

# Define arguments fill the wildcards in the below paths 
dg_first.inputs.template_args = dict(
    struct=[["subject_id","subject_id"]],
    func=[["subject_id","subject_id","run_id"]],
    events=[["subject_id","subject_id","run_id"]]
)

dg_first.inputs.field_template = dict(
    struct = "sub-%s/anat/sub-%s_T1w.nii.gz",
    func="sub-%s/func/sub-%s_task-flanker_run-%d_bold.nii.gz",
    events="sub-%s/func/sub-%s_task-flanker_run-%d_events.tsv", 
)

wf_first.connect([
        (infosource_first, dg_first, [("subject_id", "subject_id")])
])

img2float = pe.MapNode(
    interface=fsl.ImageMaths(out_data_type='float', op_string='', suffix='_dtype'),
    name='img2float',
    iterfield=['in_file'])
wf_first.connect(dg_first, 'func', img2float, 'in_file')

# Extract one roi
extract_ref = pe.MapNode(interface=fsl.ExtractROI(t_size=1), 
                         name='extractref',
                         iterfield=['in_file'])
wf_first.connect(img2float, 'out_file', extract_ref, 'in_file')

# Define a function to return the 1 based index of the middle volume
def getmiddlevolume(func):
    from nibabel import load
    funcfile = func
    if isinstance(func, list):
        funcfile = func[0]
    _, _, _, timepoints = load(funcfile).shape
    return int(timepoints/2) 

wf_first.connect(img2float, ('out_file', getmiddlevolume), extract_ref, 't_min')

motion_correct = pe.MapNode(
    interface=fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True, interpolation='spline'),
    name='realign',
    iterfield=['in_file','ref_file'])
wf_first.connect(img2float, 'out_file', motion_correct, 'in_file')
wf_first.connect(extract_ref, 'roi_file', motion_correct, 'ref_file')

plot_motion = pe.MapNode(
    interface=fsl.PlotMotionParams(in_source='fsl'),
    name='plot_motion',
    iterfield=['in_file'])
plot_motion.iterables = ('plot_type', ['rotations', 'translations','displacement'])
wf_first.connect(motion_correct, 'par_file', plot_motion, 'in_file')

meanfunc = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    name='meanfunc',
    iterfield=['in_file'])
wf_first.connect(motion_correct, 'out_file', meanfunc, 'in_file')

meanfuncmask = pe.MapNode(
    interface=fsl.BET(mask=True, no_output=True, frac=0.3),
    name='meanfuncmask',
    iterfield=['in_file'])
wf_first.connect(meanfunc, 'out_file', meanfuncmask, 'in_file')

maskfunc = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_bet', op_string='-mas'),
    name='maskfunc',
    iterfield=['in_file','in_file2'])
wf_first.connect(motion_correct, 'out_file', maskfunc, 'in_file')
wf_first.connect(meanfuncmask, 'mask_file', maskfunc, 'in_file2')

getthresh = pe.MapNode(
    interface=fsl.ImageStats(op_string='-p 2 -p 98'),
    name='getthreshold',
    iterfield=['in_file'])
wf_first.connect(maskfunc, 'out_file', getthresh, 'in_file')

threshold = pe.MapNode(
    interface=fsl.ImageMaths(out_data_type='char', suffix='_thresh'),
    name='threshold',
    iterfield=['in_file','op_string'])
wf_first.connect(maskfunc, 'out_file', threshold, 'in_file')

# Define a function to get 10% of the intensity
def getthreshop(thresh):
    return ['-thr %.10f -Tmin -bin' % (0.1 * th[1]) for th in thresh]

wf_first.connect(getthresh, ('out_stat', getthreshop), threshold, 'op_string')

medianval = pe.MapNode(
    interface=fsl.ImageStats(op_string='-k %s -p 50'),
    name='medianval',
    iterfield=['in_file','mask_file'])
wf_first.connect(motion_correct, 'out_file', medianval, 'in_file')
wf_first.connect(threshold, 'out_file', medianval, 'mask_file')

dilatemask = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_dil', op_string='-dilF'),
    name='dilatemask',
    iterfield=['in_file'])
wf_first.connect(threshold, 'out_file', dilatemask, 'in_file')

maskfunc2 = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_thresh', op_string='-mas'),
    iterfield=['in_file','in_file2'],
    name='maskfunc2')
wf_first.connect(motion_correct, 'out_file', maskfunc2, 'in_file')
wf_first.connect(dilatemask, 'out_file', maskfunc2, 'in_file2')

meanfunc2 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    name='meanfunc2',
    iterfield=['in_file'])
wf_first.connect(maskfunc2, 'out_file', meanfunc2, 'in_file')

# why here use node not mapnode?
mergenode = pe.Node(interface=util.Merge(2, axis='hstack'), 
                       name='merge')
wf_first.connect(meanfunc2, 'out_file', mergenode, 'in1')
wf_first.connect(medianval, 'out_stat', mergenode, 'in2')

fwhm_thr = 4.9996179300001655

smooth = pe.MapNode(
    interface=fsl.SUSAN(fwhm = fwhm_thr),
    name='smooth',
    iterfield=['in_file', 'brightness_threshold', 'usans'])

# Define a function to get the brightness threshold for SUSAN
def getbtthresh(medianvals):
    return [0.75 * val for val in medianvals]

def getusans(x):
    return [[tuple([val[0], 0.75 * val[1]])] for val in x]

wf_first.connect(maskfunc2, 'out_file', smooth, 'in_file')
wf_first.connect(medianval, ('out_stat', getbtthresh), smooth,
                'brightness_threshold')
wf_first.connect(mergenode, ('out', getusans), smooth, 'usans')

maskfunc3 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-mas'),
    name='maskfunc3',
    iterfield=['in_file','in_file2'])
wf_first.connect(smooth, 'smoothed_file', maskfunc3, 'in_file')
wf_first.connect(dilatemask, 'out_file', maskfunc3, 'in_file2')

intnorm = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_intnorm'),
    iterfield=['in_file', 'op_string'],
    name='intnorm')
wf_first.connect(maskfunc3, 'out_file', intnorm, 'in_file')

# Define a function to get the scaling factor for intensity normalization
def getinormscale(medianvals):
    return ['-mul %.10f' % (10000. / val) for val in medianvals]

wf_first.connect(medianval, ('out_stat', getinormscale), intnorm, 'op_string')

# Generate a mean functional image from the scaled data and this mean func will be used in performing temporal highpass filtering
meanfunc3 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean'),
    iterfield=['in_file'],
    name='meanfunc3')
wf_first.connect(intnorm, 'out_file', meanfunc3, 'in_file')

# Perform temporal highpass filtering on the data
hpcutoff = 100
TR = 2.  # ensure float

highpass = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_tempfilt'),
    name='highpass',
    iterfield=['in_file','op_string'])

# 25 = (hpcutoff / 2*TR) not (hpcutoff / TR)
def gethpstring(tempMean):
    return ['-bptf 25 -1 -add %s' % (tm) for tm in tempMean]

wf_first.connect(intnorm, 'out_file', highpass, 'in_file')
wf_first.connect(meanfunc3, ('out_file',gethpstring), highpass, 'op_string')

meanfunc4 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    iterfield=['in_file'],
    name='meanfunc4')
wf_first.connect(highpass, 'out_file', meanfunc4, 'in_file')

def subjinfo(events):
    from nipype.interfaces.base import Bunch
    import pandas as pd
    import numpy as np

    subject_info = []

    ev = pd.read_csv(events, sep="\t")
    ev = ev[ev['correctness']=='correct']
    ev['new_type'] = ev['trial_type'].apply(lambda x: str(x).split('_')[0])

    run_info = Bunch(onsets=[], 
                     durations=[])

    run_info.set(conditions=[g[0] for g in ev.groupby("new_type")])

    for group in ev.groupby("new_type"):
        run_info.onsets.append(group[1].onset.tolist())
        run_info.durations.append(group[1].duration.tolist())
    subject_info.append(run_info)

    return subject_info

get_sub_info = pe.MapNode(
    Function(
        function=subjinfo, input_names=["events"], output_names="subject_info"
    ),
    name="get_sub_info", iterfield=["events"]
)

# Connect to workflow
wf_first.connect(dg_first, 'events', get_sub_info, 'events')

condition_names = ["congruent", "incongruent"]

# Activation Contrasts (similar to https://direct.mit.edu/jocn/article/23/10/3162/5327/Is-Morality-Unified-Evidence-that-Distinct-Neural)
## From FEAT: "The correct way to tell whether two conditions are significantly different is to run a differential contrast like [1 -1] between them"
cont01 = ("congruent", "T", condition_names,  [1, 0])
cont02 = ("incongruent", "T", condition_names, [0, 1])
cont03 = ("congruent-incongruent", "T", condition_names, [1, -1])
cont04 = ("incongruent-congruent", "T", condition_names,  [-1, 1])

contrast_list = [
    cont01,
    cont02,
    cont03,
    cont04
]

TR = 2.0  # Repetition Time
hpcutoff = 100  # highpass filter cutoff in seconds! 

modelspec = pe.MapNode(interface=model.SpecifyModel(), name="modelspec", iterfield=["subject_info","functional_runs"])
modelspec.inputs.input_units = "secs"
modelspec.inputs.high_pass_filter_cutoff = hpcutoff
modelspec.inputs.time_repetition = TR

wf_first.connect(get_sub_info, 'subject_info', modelspec,'subject_info')
wf_first.connect(highpass, 'out_file',modelspec,'functional_runs')

level1design = pe.MapNode(interface=fsl.Level1Design(), name="level1design", iterfield=["session_info"])
level1design.inputs.interscan_interval = TR
# Set HRF bases functions
level1design.inputs.bases = {"gamma": {"derivs": False, 'gammasigma':3, 'gammadelay':6}}
level1design.inputs.model_serial_correlations = True
level1design.inputs.contrasts = contrast_list

wf_first.connect(modelspec,'session_info', level1design, 'session_info')

modelgen = pe.MapNode(
    interface=fsl.FEATModel(),
    name='modelgen',
    iterfield=['fsf_file', 'ev_files'])

wf_first.connect([(level1design, modelgen, [('fsf_files', 'fsf_file'), 
                                    ('ev_files','ev_files')])])

level1estimate = pe.MapNode(
    interface=fsl.FILMGLS(smooth_autocorr=True, mask_size=5, threshold=1000),
    name='level1estimate',
    iterfield=['design_file', 'in_file', 'tcon_file'])

wf_first.connect([
    (highpass,level1estimate,[('out_file','in_file')]),
    (modelgen,level1estimate,[('design_file','design_file'),
                               ('con_file','tcon_file')])
])

flt = pe.MapNode(interface=fsl.FLIRT(cost='corratio', dof=12,
                                    searchr_x = [-90,90],
                                    searchr_y = [-90,90],
                                    searchr_z = [-90,90],
                                    interp = 'trilinear'), 
               name='example_func2standard',
               iterfield=['in_file'])

flt.inputs.reference = mni_brain #fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
wf_first.connect(extract_ref, 'roi_file', flt, 'in_file')

convertxfm = pe.MapNode(interface=fsl.ConvertXFM(invert_xfm = True), 
                         name='convertxfm',
                         iterfield=['in_file'])
wf_first.connect(flt, 'out_matrix_file', convertxfm, 'in_file')

def warp_files(copes, varcopes, masks, mat, reference):
    # need to reimport here, otherwise errors come out
    import nipype.interfaces.fsl as fsl 
    
    out_copes = []
    out_varcopes = []
    out_masks = []
    
    # register mask, same function, different parameters
    warp_mask = fsl.FLIRT(apply_xfm = True, 
                     interp = 'nearestneighbour')
    warp_mask.inputs.reference = reference #fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
    warp_mask.inputs.in_matrix_file = mat
    warp_mask.inputs.output_type = "NIFTI_GZ"
    warp_mask.inputs.in_file = masks
    res_mask = warp_mask.run()
    out_masks.append(str(res_mask.outputs.out_file))
    
    # register copes & varcopes using same function, different parameters
    warp = fsl.FLIRT(apply_xfm = True, 
                     interp = 'trilinear')
    warp.inputs.reference = reference #fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
    warp.inputs.in_matrix_file = mat
    warp.inputs.output_type = "NIFTI_GZ"
    
    # register copes
    for cope in copes:
        warp.inputs.in_file = cope
        res = warp.run()
        out_copes.append(str(res.outputs.out_file))
        
     # register varcopes
    for varcope in varcopes:
        warp.inputs.in_file = varcope
        res = warp.run()
        out_varcopes.append(str(res.outputs.out_file))

    return out_copes, out_varcopes, out_masks

warpfunc = pe.MapNode(util.Function(input_names=['copes', 'varcopes', 'masks', 'mat', 'reference'],
                                    output_names=['out_copes', 'out_varcopes', 'out_masks'],
                                    function=warp_files),
                      iterfield=['copes', 'varcopes', 'masks','mat'],
                      name='warpfunc')
warpfunc.inputs.reference = mni_brain

wf_first.connect(level1estimate, 'copes', warpfunc, 'copes')
wf_first.connect(level1estimate, 'varcopes', warpfunc, 'varcopes')
wf_first.connect(dilatemask, 'out_file', warpfunc, 'masks')
wf_first.connect(flt, 'out_matrix_file', warpfunc, 'mat')

# save all results into single/separate folder for further analysis use
datasink_first = pe.Node(nio.DataSink(), name='sinker')
datasink_first.inputs.base_directory=os.path.join(exp_dir_first, "level1_results")

wf_first.connect(infosource_first, 'subject_id', datasink_first, 'container')
wf_first.connect([(level1estimate, datasink_first, [('results_dir', 'results_dir')])])
wf_first.connect([(warpfunc, datasink_first, [('out_copes', 'reg_copes')])])
wf_first.connect([(warpfunc, datasink_first, [('out_varcopes', 'reg_varcopes')])])
wf_first.connect([(warpfunc, datasink_first, [('out_masks', 'reg_masks')])])

# Create 1st-level analysis output graph
wf_first.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_first.base_dir, wf_first.name, 'graph.png'))

import logging
logging.getLogger('nipype.workflow').setLevel(0) #set to 1 to get nipype output
# Run Workflow
wf_first.run(plugin="MultiProc", plugin_args={"n_procs": 8})

fsaverage = nilearn.datasets.fetch_surf_fsaverage()
plotting.show() 

nipype_cope1 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope1.nii.gz'
nipype_cope2 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope2.nii.gz'
nipype_cope3 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope3.nii.gz'
nipype_cope4 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope4.nii.gz'

plotting.plot_stat_map(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope2, bg_img=nipype_cope1, title = 'Nipype - COPE2', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope3, bg_img=nipype_cope1, title = 'Nipype - COPE3', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope4, bg_img=nipype_cope1, title = 'Nipype - COPE4', cmap = 'bwr', colorbar = False)
plt.show()

wf_second = pe.Workflow(name='level2', base_dir=exp_dir_second)
wf_second.config["execution"]["crashfile_format"] = "txt"

# we want to group the outcome by contrast not subject
contr_list = [1,2,3,4]
infosource_second = pe.Node(util.IdentityInterface(fields=["contr_id"]),
                  name="infosource_second")
infosource_second.iterables = [("contr_id", contr_list)]

# here we use SelectFiles, instead of DataGrabber, because the former is more flexible with formatting syntax
templates_second = {
    "reg_copes":"*/reg_copes/*/*/cope{contr_id}_flirt.nii.gz",
    "reg_varcopes":"*/reg_varcopes/*/*/varcope{contr_id}_flirt.nii.gz",
    "reg_masks":"*/reg_masks/*/*/*.nii.gz"
}

dg_second = pe.Node(interface=nio.SelectFiles(templates_second),
             name="selectfiles_second")
dg_second.inputs.base_directory = data_dir_second

wf_second.connect([
        (infosource_second, dg_second, [("contr_id", "contr_id")])
])

copemerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="copemerge_second")

varcopemerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="varcopemerge_second")

maskmerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="maskmerge_second")


def sort_copes(files):
    numelements = len(files[0])
    outfiles = []
    for i in range(numelements):
        outfiles.insert(i, [])
        for j, elements in enumerate(files):
            outfiles[i].append(elements[i])
    return outfiles

wf_second.connect(dg_second, 'reg_varcopes', varcopemerge_second, 'in_files')
wf_second.connect(dg_second, 'reg_copes', copemerge_second, 'in_files')
wf_second.connect(dg_second, 'reg_masks', maskmerge_second, 'in_files')

minmask_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-Tmin'),
    iterfield=['in_file'],
    name='minmask_second')

wf_second.connect(maskmerge_second, 'merged_file', minmask_second, 'in_file')

maskcope_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskcope_second')

maskvarcope_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskvarcope_second')

wf_second.connect(copemerge_second, 'merged_file', maskcope_second, 'in_file')
wf_second.connect(minmask_second, 'out_file', maskcope_second, 'in_file2')
wf_second.connect(varcopemerge_second, 'merged_file', maskvarcope_second, 'in_file')
wf_second.connect(minmask_second, 'out_file', maskvarcope_second, 'in_file2')

def get_contrasts_l2(in_files):
    import numpy as np
    import glob
    total = len(in_files)
    print(in_files)
    print(total)
    n_sub = 9  
    ev_list = ['ev'+str(x) for x in range(1,n_sub+1)]
    weight_mtx = np.zeros((n_sub,n_sub))
    weight_mtx = weight_mtx.astype(np.float64)
    np.fill_diagonal(weight_mtx,1.)
    #contr = ['','T',ev_list,list(weight_mtx[0])] --> original paper
    contr = np.array(['','T', ev_list, list(weight_mtx[0])], dtype=object) #adaptation so the construction of the array works
    contr_lst = np.tile(contr, (n_sub,n_sub))
    contr_lst = [tuple(x) for x in contr_lst] 
    for i in range(n_sub):
        # Instead of trying to change the tuple, create a new tuple with the changed value.
        contr_lst[i] = (contr_lst[i][0], contr_lst[i][1], contr_lst[i][2], list(weight_mtx[i]))

    reg_dict = {k:None for k in ev_list}
    for k in reg_dict.keys():
        start_lst = [0.0] * total
        idx = ev_list.index(k)
        start_idx = idx*2
        end_idx = idx*2 + 1
        start_lst[start_idx] = 1.
        start_lst[end_idx] = 1.
        reg_dict[k] = start_lst    
    return contr_lst, reg_dict

contrastgen_l2 = pe.Node(util.Function(input_names=['in_files'],
                                    output_names=['contr_lst', 'reg_dict'],
                                    function=get_contrasts_l2),
                      iterfield=['in_files'],
                      name='contrastgen_l2')

wf_second.connect(dg_second,'reg_copes', contrastgen_l2, 'in_files')

level2model = pe.Node(interface=fsl.MultipleRegressDesign(),
                      name='l2model')

wf_second.connect([(contrastgen_l2, level2model, [('contr_lst','contrasts'),
                                        ('reg_dict','regressors')])])

level2estimate = pe.Node(
    interface=fsl.FLAMEO(run_mode='fe'),
    name="level2estimate",
    iterfield=['cope_file', 'var_cope_file'])

wf_second.connect([
    (maskcope_second, level2estimate, [('out_file', 'cope_file')]),
    (maskvarcope_second, level2estimate, [('out_file', 'var_cope_file')]),
    (minmask_second, level2estimate, [('out_file', 'mask_file')]),
    (level2model, level2estimate, [('design_mat', 'design_file'),
                                   ('design_con', 't_con_file'), 
                                   ('design_grp', 'cov_split_file')]),
])

datasink_second = pe.Node(nio.DataSink(), name='sinker')
datasink_second.inputs.base_directory=os.path.join(exp_dir_second, "level2_results")

int2string = lambda x: 'contrast_'+str(x)
    
wf_second.connect(infosource_second, ('contr_id',int2string), datasink_second, 'container')
wf_second.connect([(level2estimate, datasink_second, [('stats_dir', 'stats_dir')])])

# Create 1st-level analysis output graph
wf_second.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_second.base_dir, wf_second.name, 'graph.png'))

# Run Workflow
wf_second.run(plugin="MultiProc", plugin_args={"n_procs": 8})

nipype_cope1 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope1_flirt.nii.gz'
nipype_cope2 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope2_flirt.nii.gz'
nipype_cope3 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope3_flirt.nii.gz'
nipype_cope4 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope4_flirt.nii.gz'

plotting.plot_stat_map(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope2, bg_img=nipype_cope2, title = 'Nipype - COPE2', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope3, bg_img=nipype_cope3, title = 'Nipype - COPE3', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope4, bg_img=nipype_cope4, title = 'Nipype - COPE4', cmap = 'bwr', colorbar = False)
plt.show()

nipype_cope1 = './output_level2/level2_results/contrast_1/stats_dir/_contr_id_1/stats/cope1.nii.gz'
nipype_cope2 = './output_level2/level2_results/contrast_2/stats_dir/_contr_id_2/stats/cope1.nii.gz'
nipype_cope3 = './output_level2/level2_results/contrast_3/stats_dir/_contr_id_3/stats/cope1.nii.gz'
nipype_cope4 = './output_level2/level2_results/contrast_4/stats_dir/_contr_id_4/stats/cope1.nii.gz'

plotting.plot_glass_brain(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1')
plotting.plot_glass_brain(nipype_cope2, bg_img=nipype_cope2, title = 'Nipype - COPE2')
plotting.plot_glass_brain(nipype_cope3, bg_img=nipype_cope3, title = 'Nipype - COPE3')
plotting.plot_glass_brain(nipype_cope4, bg_img=nipype_cope4, title = 'Nipype - COPE4')
plt.show()

wf_third = pe.Workflow(name='level3', base_dir=exp_dir_third)
wf_third.config["execution"]["crashfile_format"] = "txt"

# we want to group the outcome by contrast not subject
contr_list = [1,2,3,4]
infosource_third = pe.Node(util.IdentityInterface(fields=["contr_id"]),
                  name="infosource_third")
infosource_third.iterables = [("contr_id", contr_list)]

# here we use SelectFiles, instead of DataGrabber, because the former is more flexible with formatting syntax
templates_third = {
    "copes":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/cope*.nii.gz",
    "varcopes":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/varcope*.nii.gz",
    "masks":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/mask.nii.gz"
}
templates = {
    "reg_copes":"*/reg_copes/*/*/cope{contr_id}_flirt.nii.gz",
    "reg_varcopes":"*/reg_varcopes/*/*/varcope{contr_id}_flirt.nii.gz",
    "reg_masks":"*/reg_masks/*/*/*.nii.gz"
}

dg_third = pe.Node(interface=nio.SelectFiles(templates_third),
             name="selectfiles_third")
dg_third.inputs.base_directory = data_dir_third

wf_third.connect([
        (infosource_third, dg_third, [("contr_id", "contr_id")])
])

copemerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="copemerge_third")

varcopemerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="varcopemerge_third")

maskmerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="maskmerge_third")

def repeat_mask(file):
    import numpy as np
    import glob
    # n_sub = len(glob.glob(os.path.join(os.getcwd(), 'ds000102')+"/sub-0*"))    
    n_sub = 9
    mask_lst = [file]
    repeated = np.repeat(mask_lst,n_sub)
    return list(repeated)

wf_third.connect(dg_third, 'copes', copemerge_third, 'in_files')
wf_third.connect(dg_third, 'varcopes', varcopemerge_third, 'in_files')
wf_third.connect(dg_third, ('masks', repeat_mask), maskmerge_third, 'in_files')

# /usr/local/fsl/bin/fslmaths mask -Tmin mask
minmask_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-Tmin'),
    iterfield=['in_file'],
    name='minmask_third')

wf_third.connect(maskmerge_third, 'merged_file', minmask_third, 'in_file')

maskcope_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskcope_third')

maskvarcope_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskvarcope_third')

wf_third.connect(copemerge_third, 'merged_file', maskcope_third, 'in_file')
wf_third.connect(minmask_third, 'out_file', maskcope_third, 'in_file2')
wf_third.connect(varcopemerge_third, 'merged_file', maskvarcope_third, 'in_file')
wf_third.connect(minmask_third, 'out_file', maskvarcope_third, 'in_file2')

def num_copes(files):
    return len(files)

level3model = pe.Node(interface=fsl.L2Model(), name='l3model')
wf_third.connect(dg_third, ('copes',num_copes), level3model, 'num_copes')

level3estimate = pe.Node(
    interface=fsl.FLAMEO(run_mode='flame1'),
    name="level3estimate",
    iterfield=['cope_file', 'var_cope_file'])

wf_third.connect([
    (maskcope_third, level3estimate, [('out_file', 'cope_file')]),
    (maskvarcope_third, level3estimate, [('out_file', 'var_cope_file')]),
    (minmask_third, level3estimate, [('out_file', 'mask_file')]),
    (level3model, level3estimate, [('design_mat', 'design_file'),
                                   ('design_con', 't_con_file'), 
                                   ('design_grp', 'cov_split_file')]),
])

smoothest = pe.Node(
    interface=fsl.SmoothEstimate(dof=25),
    name="smoothest",
    iterfield=['residual_fit_file', 'mask_file'])

wf_third.connect(minmask_third, 'out_file', smoothest, 'mask_file')
wf_third.connect(level3estimate, 'res4d', smoothest, 'residual_fit_file')

level3mask = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file','in_file2'],
    name='level3mask')

wf_third.connect([
    (level3estimate,level3mask,[('zstats', 'in_file')]),
    (minmask_third, level3mask, [('out_file', 'in_file2')]),
])

level3threshold = pe.Node(
        interface=fsl.Cluster(threshold = 3.1,
                              pthreshold = 0.05,
                              use_mm = True,
                              out_threshold_file =True,
                              out_index_file = True,
                              out_localmax_txt_file = True),
    iterfield=['in_file','cope_file'],
    name='level3threshold')



wf_third.connect([
    (level3mask, level3threshold, [('out_file', 'in_file')]),
    (level3estimate, level3threshold, [('copes', 'cope_file')]),
    (smoothest, level3threshold, [('dlh', 'dlh'),
                                   ('volume', 'volume')]),
])

datasink_third = pe.Node(nio.DataSink(), name='sinker')
datasink_third.inputs.base_directory=os.path.join(exp_dir_third, "level3_results")

int2string = lambda x: 'contrast_'+str(x)
    
wf_third.connect(infosource_third, ('contr_id',int2string), datasink_third, 'container')
wf_third.connect([
    (level3estimate, datasink_third, [('stats_dir', 'stats_dir')]),
    (level3threshold, datasink_third, [('threshold_file', 'thresholded'),
                                ('localmax_txt_file', 'localmax_txt'),
                                ('index_file', 'index')]),
           ])

# Create 1st-level analysis output graph
wf_third.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_third.base_dir, wf_third.name, 'graph.png'))

# Run Workflow
logging.getLogger('nipype.workflow').setLevel(0) 
wf_third.run(plugin="MultiProc", plugin_args={"n_procs": 8})

nipype_cope1 = './output_level3/level3_results/contrast_1/stats_dir/_contr_id_1/stats/cope1.nii.gz'
nipype_cope2 = './output_level3/level3_results/contrast_2/stats_dir/_contr_id_2/stats/cope1.nii.gz'
nipype_cope3 = './output_level3/level3_results/contrast_3/stats_dir/_contr_id_3/stats/cope1.nii.gz'
nipype_cope4 = './output_level3/level3_results/contrast_4/stats_dir/_contr_id_4/stats/cope1.nii.gz'

t_nipype_cope1 = surface.vol_to_surf(nipype_cope1, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope1, hemi='right',
        title="Nipype - COPE1", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope2 = surface.vol_to_surf(nipype_cope2, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope2, hemi='right',
        title="Nipype - COPE2", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope3 = surface.vol_to_surf(nipype_cope3, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope3, hemi='right',
        title="Nipype - COPE3", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope4 = surface.vol_to_surf(nipype_cope4, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope4, hemi='right',
        title="Nipype - COPE4", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()

nipype_cope1 = './output_level3/level3_results/contrast_1/thresholded/_contr_id_1/zstat1_maths_threshold.nii.gz'
nipype_cope2 = './output_level3/level3_results/contrast_2/thresholded/_contr_id_2/zstat1_maths_threshold.nii.gz'
nipype_cope3 = './output_level3/level3_results/contrast_3/thresholded/_contr_id_3/zstat1_maths_threshold.nii.gz'
nipype_cope4 = './output_level3/level3_results/contrast_4/thresholded/_contr_id_4/zstat1_maths_threshold.nii.gz'

plotting.plot_img_on_surf(nipype_cope1,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats1",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope2,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats2",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope3,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats3",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope4,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats4",
                          cmap = 'hot', colorbar=True)
plotting.show()

import os

%load_ext watermark

%watermark
%watermark --iversions

neurodesktop_version = (
    os.environ.get('JUPYTER_IMAGE', '').split(':')[-1] or
    os.environ.get('NEURODESKTOP_VERSION', 'unknown')
)

print(f"Neurodesktop version: {neurodesktop_version}")