<jupyter_start><jupyter_text>1st-level Analysis    This notebook performes the 1st-level analysis in subject space by executing the following steps:1. Aggregate 1st-level model parameters2. Specify 1st-level contrasts3. Estimate 1st-level contrasts4. Normalization to template space (optional)The notebook runs on all functional runs of a particular task, computes the specified beta-contrast and normalizes them into template space. If requested, it will also compute one beta contrast per condition, per run (usefull for a multivariate approache). Confounding factors and outlier volumes can be added as nuisance regressors in the GLM.**Note:** This notebook requires that the functional preprocessing pipeline was already executed and that it's output can be found in the dataset folder under `dataset/derivatives/fmriflows/preproc_func`.   Data Structure RequirementsThe data structure to run this notebook should be according to the BIDS format:    dataset    ├── fmriflows_spec_analysis.json    ├── sub-{sub_id}    │   └── func    │       └── sub-{sub_id}_task-{task_id}[_run-{run_id}]_events.tsv    └── derivatives        └── fmriflows            ├── preproc_anat            │   └── sub-{sub_id}            │       └── sub-{sub_id}_transformComposite.h5            └── preproc_func                └── sub-{sub_id}                    ├── sub-{sub_id}_task-{task_id}[_run-{run_id}]_nss.txt                    ├── sub-{sub_id}_task-{task_id}[_run-{run_id}]_tFilter_*_confounds.tsv                    ├── sub-{sub_id}_task-{task_id}[_run-{run_id}]_tFilter_*_confounds_outliers.txt                    └── sub-{sub_id}_task-{task_id}[_run-{run_id}]_tFilter_*_sFilter_*.nii.gz                    **Note:** `Session` and `run` identifiers are optional.`fmriflows` will run the 1st-level analysis on all runs, but separately for each session and particular task.  Execution SpecificationsThis notebook will extract the relevant analysis specifications from the `fmriflows_spec_analysis.json` file in the dataset folder. In the current setup, they are as follows:<jupyter_code>import json
from os.path import join as opj

spec_file = opj('/data', 'fmriflows_spec_analysis.json')

with open(spec_file) as f:
    specs = json.load(f)
# Extract parameters for 1st-level analysis workflow
subject_list = specs['subject_list']
session_list = specs['session_list']
tasks = specs['tasks']
filters_spatial = specs['filters_spatial']
filters_temporal = specs['filters_temporal']
nuisance_regressors = specs['nuisance_regressors']
use_outliers = specs['use_outliers']
model_serial_correlations = specs['model_serial_correlations']
model_bases = specs['model_bases']
estimation_method = specs['estimation_method']
normalize = specs['normalize']
norm_res = specs['norm_res']
con_per_run = specs['con_per_run']
norm_res_multi = specs['norm_res_multi']
postfix = specs['analysis_postfix']
n_proc = specs['n_parallel_jobs']<jupyter_output><empty_output><jupyter_text>If you'd like to change any of those values manually, overwrite them below:<jupyter_code># List of subject identifiers
subject_list
# List of session identifiers
session_list
# List of spatial filters (smoothing) that were used during functional preprocessing
filters_spatial
# List of temporal filters that were used during functional preprocessing
filters_temporal
# Nuisance identifiers that should be included in the GLM
nuisance_regressors
# If outliers detected during functional preprocing should be used in GLM
use_outliers
# Serial Correlation Model to use: 'AR(1)', 'FAST' or 'none'
model_serial_correlations
# Model bases to use: 'hrf', 'fourier', 'fourier', 'fourier_han', 'gamma' or 'fir'
# If 'hrf', also specify if time and dispersion derivatives should be used
model_bases
# Estimation Method to use: 'Classical', 'Bayesian' or 'Bayesian2'
estimation_method
# Specify if contrasts should be normalized to template space after estimation
normalize
# Specify voxel resolution of normalized contrasts
norm_res
# Specify if a contrast should be computed for stimuli category per run
con_per_run
# Specify voxel resolution of normalized contrasts for multivariate analysis
norm_res_multi
# Specify a particular analysis postfix
postfix
# Number of parallel jobs to run
n_proc
# Task specific parameters
tasks
norm_func = False<jupyter_output><empty_output><jupyter_text>Create the WorkflowTo ensure a good overview of the 1st-level analysis, the workflow was divided into an analysis and a report subworkflow. Import Modules<jupyter_code>from os.path import join as opj
from nipype import Node, MapNode, Workflow
from nipype.interfaces.utility import Function, IdentityInterface
from nipype.algorithms.misc import Gunzip
from nipype.algorithms.modelgen import SpecifySPMModel
from nipype.interfaces.spm import Level1Design, EstimateModel, EstimateContrast
from nipype.interfaces.ants import ApplyTransforms
from nipype.interfaces.utility import Merge
from nipype.interfaces.io import SelectFiles, DataSink
# Specify SPM location
from nipype.interfaces.matlab import MatlabCommand
MatlabCommand.set_default_paths('/opt/spm12-r7219/spm12_mcr/spm12')<jupyter_output><empty_output><jupyter_text>Relevant Execution Variables<jupyter_code># Folder paths and names
exp_dir = '/data/derivatives'
out_dir = 'fmriflows'
work_dir = '/workingdir'<jupyter_output><empty_output><jupyter_text>Create a subworkflow for the Analysis Workflow<jupyter_code># Specify 1st-level model parameters (stimuli onsets, duration, etc.)
def collect_model_info(event_files, nss_files, condition_names, TR):

    import numpy as np
    import pandas as pd
    from nipype.interfaces.base import Bunch

    model_info = []
    stimuli_order = []

    for i, f in enumerate(event_files):

        trialinfo = pd.read_csv(f, sep='	')
        stimuli_list = [t for t in trialinfo.trial_type if str(t) in condition_names]
        stimuli_order.append(stimuli_list)
        nss = np.loadtxt(nss_files[i])
        conditions = []
        onsets = []
        durations = []

        for group in trialinfo.groupby('trial_type'):
            if str(group[0]) in condition_names:
                conditions.append(str(group[0]))
                onsets.append(list(group[1].onset - nss * TR))
                durations.append(group[1].duration.tolist())

        model_info.append(Bunch(conditions=conditions,
                                onsets=onsets,
                                durations=durations))
   
    return model_info, stimuli_order

# Get Subject Info - get subject specific condition information
get_model_info = Node(Function(input_names=['event_files', 'nss_files', 'condition_names', 'TR'],
                               output_names=['model_info', 'stimuli_order'],
                               function=collect_model_info),
                      name='get_model_info')
# Gunzip NIfTI files for SPM
gunzip = MapNode(Gunzip(), name='gunzip', iterfield=['in_file'])
# Create nuisance regressors
def create_nuisance_regressors(confounds, nuisance_regressors):

    import numpy as np
    import pandas as pd
    from os.path import basename, abspath

    # To store regressor files into
    regressor_files = []

    # Go through confound files
    for i, c in enumerate(confounds):
        df = pd.read_csv(c, sep='	')
        selection = [k for k in df.keys() for n in nuisance_regressors if n in k]
        dfs = df[selection]
        out_file = abspath(basename('confounds_%02d.rst' % (i + 1)))
        np.savetxt(out_file, dfs.values)
        regressor_files.append(out_file)

    return regressor_files

nuisance_reg = Node(Function(input_names=['confounds', 'nuisance_regressors'],
                             output_names=['confounds'],
                             function=create_nuisance_regressors),
                      name='nuisance_reg')
nuisance_reg.inputs.nuisance_regressors = nuisance_regressors
# Create SPM model
model_spec = Node(SpecifySPMModel(concatenate_runs=False,
                                  input_units='secs',
                                  output_units='secs'),
                  name="model_spec")
# Create 1st-level desing
level1_design = Node(Level1Design(bases=model_bases,
                                  timing_units='secs',
                                  model_serial_correlations=model_serial_correlations),
                    name="level1_design")
# Estimate 1st-level model
level1_estimate = Node(EstimateModel(estimation_method=estimation_method),
                       name="level1_estimate")
# Estimate 1st-level contrasts
level1_con_est = Node(EstimateContrast(), name="level1_con_est")
# Create analysis workflow
analysisflow = Workflow(name='analysisflow')

# Add nodes to workflow and connect them
analysisflow.connect([(get_model_info, model_spec, [('model_info', 'subject_info')]),
                      (gunzip, model_spec, [('out_file', 'functional_runs')]),
                      (nuisance_reg, model_spec, [('confounds', 'realignment_parameters')]),
                      (model_spec, level1_design, [('session_info', 'session_info')]),
                      (level1_design, level1_estimate, [('spm_mat_file', 'spm_mat_file')]),
                      (level1_estimate, level1_con_est, [('spm_mat_file', 'spm_mat_file'),
                                                         ('beta_images', 'beta_images'),
                                                         ('residual_image', 'residual_image')]),
                     ])<jupyter_output><empty_output><jupyter_text>Create a subworkflow for the Report Workflow<jupyter_code># Plots design matrix
def plot_design_matrix(SPM):

    import numpy as np
    from matplotlib import pyplot as plt
    from scipy.io import loadmat
    from os.path import basename, abspath

    # Using scipy's loadmat function we can access SPM.mat
    spmmat = loadmat(SPM, struct_as_record=False)
    
    # Now we can load the design matrix and the names of the rows
    designMatrix = spmmat['SPM'][0][0].xX[0][0].X
    names = [i[0] for i in spmmat['SPM'][0][0].xX[0][0].name[0]]

    # Value normalization for better visualization
    normed_design = designMatrix / np.abs(designMatrix).max(axis=0)

    # Plotting of the design matrix
    fig, ax = plt.subplots(figsize=(8, 8))
    plt.imshow(normed_design, aspect='auto', cmap='gray', interpolation='nearest')
    ax.set_ylabel('Volume id')
    ax.set_xticks(np.arange(len(names)))
    ax.set_xticklabels(names, rotation=90)
    design_matrix = abspath(basename(SPM.replace('.mat', '.png')))
    
    fig.tight_layout()
    fig.savefig(design_matrix)
    
    return design_matrix

# Extracts design matrix from SPM.mat and plots it
plot_GLM = Node(Function(input_names=['SPM'],
                         output_names=['out_file'],
                         function=plot_design_matrix),
                name='plot_GLM')
# Plot contrast image
def plot_contrast(con_files, contrast_list, threshold):

    from nilearn.image import math_img
    from nilearn.plotting import plot_glass_brain
    from os.path import basename, abspath

    out_files = []
    for i, c in enumerate(con_files):
        percentile_str = 'np.percentile(np.abs(img[np.nan_to_num(img)!=0]), %s)' % threshold
        con_thr = math_img('img * (np.abs(np.nan_to_num(img))>=%s)' % percentile_str, img=c)
        title = contrast_list[i][0]
        new_file_name = c.replace('.nii', '.png').replace('.gz', '')
        out_file = abspath(basename(new_file_name))
        plot_glass_brain(con_thr, display_mode='lyrz', colorbar=True,
                         symmetric_cbar=False, plot_abs=False, black_bg=True,
                         title=title, output_file=out_file)
        out_files.append(out_file)
    
    return out_files

# Extracts design matrix from SPM.mat and plots it
plot_contrasts = Node(Function(input_names=['con_files', 'contrast_list', 'threshold'],
                               output_names=['out_files'],
                               function=plot_contrast),
                      name='plot_contrasts')
plot_contrasts.inputs.threshold = 98
# Create report workflow
reportflow = Workflow(name='reportflow')

# Add nodes to workflow and connect them
reportflow.add_nodes([plot_GLM,
                      plot_contrasts,
                     ])<jupyter_output><empty_output><jupyter_text>Specify Input & Output Stream<jupyter_code># Iterate over subject, session, task and run id
info_source = Node(IdentityInterface(fields=['subject_id',
                                             'session_id',
                                             'task_id',
                                             'spatial_filt',
                                             'temporal_filt']),
                   name='info_source')

iter_list = [('subject_id', subject_list),
             ('task_id', list(tasks.keys())),
             ('spatial_filt', filters_spatial),
             ('temporal_filt', filters_temporal),
             ]

if session_list:
    iter_list.append(('session_id', session_list))
else:
    info_source.inputs.session_id = ''

info_source.iterables = iter_list
# Create path to input files
def create_file_path(subject_id, session_id, task_id, tFilter, sFilter):

    # Collect normalization matrix
    template_anat = '/data/derivatives/fmriflows/preproc_anat/sub-{0}/sub-{0}_'
    if session_id:
        template_anat += 'ses-%s_' % session_id
    template_anat += '{1}.h5'
    transforms = template_anat.format(subject_id, 'transformComposite')
    
    # Collect outputs from the functional preprocessing
    template_func = '/data/derivatives/fmriflows/preproc_func/sub-{0}/sub-{0}_'
    if session_id:
        template_func += 'ses-%s_' % session_id
    template_func += 'task-{1}*_{2}.{3}'
    
    tFilter_id = 'tFilter_%s.%s' % (tFilter[0], tFilter[1])
    sFilter_id = 'sFilter_%s_%s' % (sFilter[0], sFilter[1])

    from glob import glob
    nss = glob(template_func.format(
        subject_id, task_id, '*nss', 'txt'))
    confounds = glob(template_func.format(
        subject_id, task_id, '%s_*confounds' % tFilter_id, 'tsv'))
    outliers = glob(template_func.format(
        subject_id, task_id, '%s_*confounds_outliers' % tFilter_id, 'txt'))
    func = glob(template_func.format(
        subject_id, task_id, '%s_*%smm' % (tFilter_id, sFilter_id), 'nii.gz'))

    # Get all anatomical files
    from bids.layout import BIDSLayout
    layout = BIDSLayout('/data/')

    # Collect event files
    search_parameters = {'task': task_id,
                         'return_type': 'file',
                         'suffix': 'events'
                        }

    if session_id:
        search_parameters['session'] = session_id

    event_files = layout.get(**search_parameters)
    if len(event_files) == 1:
        events = [event_files[0]] * len(func)
    elif len(event_files) == 0:
        search_parameters['session'] = None
        event_files = layout.get(**search_parameters)
        events = [event_files[0]] * len(func)
    else:
        events = []
        for e in event_files:
            if 'sub-%s' % subject_id in e:
                events.append(e)

    return transforms, sorted(func), sorted(events), sorted(outliers), sorted(confounds), sorted(nss)

select_files = Node(Function(input_names=['subject_id', 'session_id', 'task_id',
                                          'tFilter', 'sFilter'],
                             output_names=['transforms', 'func', 'events',
                                           'outliers', 'confounds', 'nss'],
                             function=create_file_path),
                    name='select_files')
# Extract sequence specifications of functional images
def get_parameters(func, tFilter, tasks, task_id):
    
    from nibabel import load
    
    # Extract TR from first functional image
    TR = load(func[0]).header.get_zooms()[3]
    
    # Specify high pass filter
    high_pass = tFilter[1] if not None else 128.

    # Extract contrasts from specification file
    contrast_list = tasks[task_id]
    contrasts = [[c[0], c[2], contrast_list['condition_names'], c[1]]
                 for c in contrast_list['contrasts']]

    return TR, high_pass, contrasts, contrast_list['condition_names']

get_param = Node(Function(input_names=['func', 'tFilter', 'tasks', 'task_id'],
                          output_names=['TR', 'high_pass', 'contrasts', 'condition_names'],
                          function=get_parameters),
                 name='get_param')
get_param.inputs.tasks = tasks
# Save relevant outputs in a datasink
datasink = Node(DataSink(base_directory=exp_dir,
                         container=out_dir),
                name='datasink')
# Apply the following naming substitutions for the datasink
if session_list:

    folder_old = ['_session_id_%s_spatial_filt_%s_subject_id_%s_task_id_%s_temporal_filt_%s/' % (
        ses, '.'.join([str(f) for f in sFilter]), sub, task, '.'.join([str(t) for t in tFilter]))
                  for sub in subject_list
                  for ses in session_list
                  for task in list(tasks.keys())
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]

    folder_new = ['sub-%s/task-%s/ses-%s/tFilter_%s_sFilter_%s/' % (
         sub, task, ses, '.'.join([str(t) for t in tFilter]), '.'.join([str(f) for f in sFilter]))
                  for sub in subject_list
                  for ses in session_list
                  for task in list(tasks.keys())
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]
else:
    
    folder_old = ['_spatial_filt_%s_subject_id_%s_task_id_%s_temporal_filt_%s/' % (
        '.'.join([str(f) for f in sFilter]), sub, task, '.'.join([str(t) for t in tFilter]))
                  for sub in subject_list
                  for task in list(tasks.keys())
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]

    folder_new = ['sub-%s/task-%s/tFilter_%s_sFilter_%s/' % (
         sub, task, '.'.join([str(t) for t in tFilter]), '.'.join([str(f) for f in sFilter]))
                  for sub in subject_list
                  for task in list(tasks.keys())
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]
    
substitutions = [z for z in zip(folder_old, folder_new)]
substitutions += [('_gzip_con%d/' % i, '') for i in range(200)]
substitutions += [('_gzip_con_run%d/' % i, '') for i in range(200)]
datasink.inputs.substitutions = substitutions<jupyter_output><empty_output><jupyter_text>Create 1st-Level Analysis Workflow<jupyter_code># Create 1st level analysis workflow
analysis_1st = Workflow(name='analysis_1st')
analysis_1st.base_dir = work_dir
out_folder = 'analysis_1stLevel'
if postfix:
    out_folder += '_%s' % postfix

# Add nodes to workflow and connect them
analysis_1st.connect([(info_source, select_files, [('subject_id', 'subject_id'),
                                                   ('session_id', 'session_id'),
                                                   ('task_id', 'task_id'),
                                                   ('spatial_filt', 'sFilter'),
                                                   ('temporal_filt', 'tFilter')]),
                      (info_source, get_param, [('task_id', 'task_id'),
                                                ('temporal_filt', 'tFilter')]),
                      (select_files, get_param, [('func', 'func')]),
                      (get_param, analysisflow, [('TR', 'model_spec.time_repetition'),
                                                 ('high_pass', 'model_spec.high_pass_filter_cutoff'),
                                                 ('TR', 'level1_design.interscan_interval'),
                                                 ('contrasts', 'level1_con_est.contrasts'),
                                                 ('condition_names', 'get_model_info.condition_names'),
                                                 ('TR', 'get_model_info.TR'),
                                                ]),
                      (get_param, reportflow, [('contrasts', 'plot_contrasts.contrast_list')]),
                      (select_files, analysisflow, [('func', 'gunzip.in_file'),
                                                    ('events', 'get_model_info.event_files'),
                                                    ('nss', 'get_model_info.nss_files'),
                                                    ('confounds', 'nuisance_reg.confounds')]),
                      
                      # Connect analysis and report workflow
                      (analysisflow, reportflow, [('level1_con_est.spm_mat_file', 'plot_GLM.SPM')]),

                      # Store analysis results in datasink
                      (analysisflow, datasink, [
                          ('level1_estimate.beta_images', '%s.univariate.@betas' % out_folder),
                          ('level1_estimate.residual_image', '%s.univariate.@ResMS' % out_folder),
                          ('level1_con_est.spm_mat_file', '%s.univariate.@spm_mat' % out_folder),
                          ('level1_con_est.con_images', '%s.univariate.@con' % out_folder),
                          ('level1_con_est.ess_images', '%s.univariate.@ess' % out_folder),
                          ('level1_con_est.spmT_images', '%s.univariate.@spmT' % out_folder),
                          ('level1_con_est.spmF_images', '%s.univariate.@spmF' % out_folder)]),

                      # Store report results in datasink
                      (reportflow, datasink, [
                          ('plot_GLM.out_file', '%s.univariate.@spm_mat_svg' % out_folder)]),
                      ])
# Add outlier parameters if requested by user
if use_outliers:
    analysis_1st.connect([(select_files, analysisflow, [('outliers', 'model_spec.outlier_files')])])<jupyter_output><empty_output><jupyter_text>Add normalization subworkflow if requested<jupyter_code># Merge node
merge = Node(Merge(2, ravel_inputs=True), name='merge')
# Creation of template brain with desired voxel resolution
template_dir = '/templates/mni_icbm152_nlin_asym_09c/'
brain_template = opj(template_dir, '1.0mm_brain.nii.gz')

# Resample template brain to desired resolution
from nibabel import load, Nifti1Image
from nilearn.image import resample_img
from nibabel.spaces import vox2out_vox

img = load(brain_template)
target_shape, target_affine = vox2out_vox(img, voxel_sizes=norm_res)
img_resample = resample_img(img, target_affine, target_shape, clip=True)
norm_template = opj(template_dir, 'template_brain_%s.nii.gz' %'_'.join([str(n) for n in norm_res]))
img_resample.to_filename(norm_template)

# Normalize contrasts if requested
norm_con = MapNode(ApplyTransforms(reference_image=norm_template,
                                   input_image_type=3,
                                   float=True,
                                   interpolation='BSpline',
                                   invert_transform_flags=[False],
                                   out_postfix='_norm'),
                   name='norm_con', iterfield=['input_image'])
# Gzip normalized contrasts
def gzip_nifti(contrast):
    
    import nibabel as nb
    from os.path import basename, abspath
    
    # Change the compression level of the NIfTI image
    nb.openers.Opener.default_compresslevel = 6

    # Save compressed contrast
    out_file = abspath(basename(contrast.replace('.nii', '.nii.gz')))
    nb.load(contrast).to_filename(out_file)

    return out_file

gzip_con = MapNode(Function(input_names=['contrast'],
                            output_names=['out_file'],
                            function=gzip_nifti),
                   name='gzip_con', iterfield=['contrast'])
# Normalize contrasts to template space if requested by user
if normalize:
    
    # Create normalization workflow
    normflow = Workflow(name='normflow')
    normflow.base_dir = work_dir
    
    # Connect nodes within normalization workflow
    normflow.connect([(merge, norm_con, [('out', 'input_image')]),
                      (norm_con, gzip_con, [('output_image', 'contrast')]),
                     ])
    
    # Connect analysis workflow to norm subworkflow
    analysis_1st.connect([(select_files, normflow, [('transforms', 'norm_con.transforms')]),
                          (analysisflow, normflow, [('level1_con_est.con_images', 'merge.in1'),
                                                    ('level1_con_est.ess_images', 'merge.in2')]),
                          (normflow, datasink, [
                              ('gzip_con.out_file', '%s.univariate.@norm_files' % out_folder)]),
                          (normflow, reportflow, [('gzip_con.out_file', 'plot_contrasts.con_files')]),
                          (reportflow, datasink, [
                              ('plot_contrasts.out_files', '%s.univariate.@con_plots' % out_folder)]),
                         ])
else:
    # Connect contrast plot node to analysis workflow
    analysis_1st.connect([(analysisflow, merge, [('level1_con_est.con_images', 'in1'),
                                                 ('level1_con_est.ess_images', 'in2')]),
                          (merge, reportflow, [('out', 'plot_contrasts.con_files')]),
                          (reportflow, datasink, [
                              ('plot_contrasts.out_files', '%s.univariate.@con_plots' % out_folder),]),
                         ])<jupyter_output><empty_output><jupyter_text>Add subworkflow for to create contrasts for multivariate analysis<jupyter_code># Create normalization workflow
multiflow = Workflow(name='multiflow')
multiflow.base_dir = work_dir
# Create contrast list for condition per run
def get_con_per_run(stimuli_order, tasks, task_id):

    import numpy as np

    # Extract condition names
    condition_names = tasks[task_id]['condition_names']

    # Aggregate event information
    event_list = []
    for i, l in enumerate(stimuli_order):
        event_list.append([z for z in zip(np.full(len(l), i), l)])
    event_info = np.reshape(event_list, (-1, 2))
    unique_contrasts = np.unique(event_info[:,1])
    n_runs = np.unique(event_info[:,0])

    # Create list of contrasts for each condition per run
    contrast_list_run = []
    n_contrasts = len(condition_names)
    n_conditions = len(n_runs)
    condition_labels = []

    for j in range(n_conditions):
        for i in range(n_contrasts):
            name = 'cont_%05d' % (1 + i + j * n_conditions)
            con_id = np.zeros(n_contrasts).tolist()
            run_id = np.zeros(n_conditions).tolist()
            con_id[i] = 1
            run_id[j] = 1
            contrast_list_run.append([
                name, 'T', condition_names, con_id, run_id])
            condition_labels.append(condition_names[i])
    
    return contrast_list_run, condition_labels

# Extracts design matrix from SPM.mat and plots it
comp_con_per_run = Node(Function(input_names=['stimuli_order', 'tasks', 'task_id'],
                                 output_names=['run_contrasts', 'condition_labels'],
                                 function=get_con_per_run),
                        name='comp_con_per_run')
comp_con_per_run.inputs.tasks = tasks
# Estimate 1st-level contrasts - one for each session
level1_con_est_run = Node(EstimateContrast(), name="level1_con_est_run")
# Resample template brain for multivariate analysis to desired resolution
img = load(brain_template)
target_shape, target_affine = vox2out_vox(img, voxel_sizes=norm_res_multi)
img_resample = resample_img(img, target_affine, target_shape, clip=True)
norm_template_multi = opj(template_dir, 'template_brain_%s.nii.gz' %'_'.join([str(n) for n in norm_res_multi]))
img_resample.to_filename(norm_template_multi)

# Normalize contrasts
norm_con_run = MapNode(ApplyTransforms(reference_image=norm_template_multi,
                                       input_image_type=3,
                                       float=True,
                                       interpolation='BSpline',
                                       invert_transform_flags=[False],
                                       out_postfix='_norm'),
                       name='norm_con_run', iterfield=['input_image'])
# Gzip normalized contrasts
gzip_con_run = MapNode(Function(input_names=['contrast'],
                                output_names=['out_file'],
                                function=gzip_nifti),
                       name='gzip_con_run', iterfield=['contrast'])
# Write label file
def write_labels_file(condition_labels, spm_mat_file):

    import numpy as np
    from os.path import basename, abspath
    label_file = abspath(basename(spm_mat_file.replace('SPM.mat',
                                                       'labels.csv')))
    np.savetxt(label_file, condition_labels, fmt='%s')

    return label_file

write_labels_run = Node(Function(input_names=['condition_labels', 'spm_mat_file'],
                                 output_names=['labels_file'],
                                 function=write_labels_file),
                        name='write_labels_run')
# Create workflow if contrasts per condition per run should be computed
if con_per_run:

    # Connect nodes within multivariate workflow
    multiflow.connect([(comp_con_per_run, level1_con_est_run, [('run_contrasts', 'contrasts')]),
                       (level1_con_est_run, norm_con_run, [('con_images', 'input_image')]),
                       (norm_con_run, gzip_con_run, [('output_image', 'contrast')]),
                       (comp_con_per_run, write_labels_run, [('condition_labels', 'condition_labels')]),
                      ])

    # Connect all nodes in this part of the workflow
    analysis_1st.connect([(info_source, multiflow, [('task_id', 'comp_con_per_run.task_id')]),
                          (select_files, multiflow, [('transforms', 'norm_con_run.transforms')]),
                          (analysisflow, multiflow, [
                              ('get_model_info.stimuli_order', 'comp_con_per_run.stimuli_order'),
                              ('level1_estimate.spm_mat_file', 'level1_con_est_run.spm_mat_file'),
                              ('level1_estimate.beta_images', 'level1_con_est_run.beta_images'),
                              ('level1_estimate.residual_image', 'level1_con_est_run.residual_image'),
                              ('level1_estimate.spm_mat_file', 'write_labels_run.spm_mat_file'),
                          ]),
                          (multiflow, datasink, [
                              ('gzip_con_run.out_file', '%s.multivariate.@norm_files_run' % out_folder),
                              ('write_labels_run.labels_file', '%s.multivariate.@norm_labels' % out_folder),
                          ]),
                         ])<jupyter_output><empty_output><jupyter_text>Visualize Workflow<jupyter_code># Create analysis_1st output graph
analysis_1st.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph in the notebook
from IPython.display import Image
Image(filename=opj(analysis_1st.base_dir, 'analysis_1st', 'graph.png'))<jupyter_output><empty_output><jupyter_text>Run Workflow<jupyter_code># Run the workflow in parallel mode
res = analysis_1st.run(plugin='MultiProc', plugin_args={'n_procs' : n_proc})
# Save workflow graph visualizations in datasink
analysis_1st.write_graph(graph2use='flat', format='png', simple_form=True)
analysis_1st.write_graph(graph2use='colored', format='png', simple_form=True)

from shutil import copyfile
copyfile(opj(analysis_1st.base_dir, 'analysis_1st', 'graph.png'),
         opj(exp_dir, out_dir,  out_folder, 'graph.png'))
copyfile(opj(analysis_1st.base_dir, 'analysis_1st', 'graph_detailed.png'),
         opj(exp_dir, out_dir, out_folder, 'graph_detailed.png'));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10577043827703127, 'fraction_numerical': 0.005619741303528821, 'mean_word_length': 2.570964125560538, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 68, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '533822', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 9718, 'n_tokens_neox': 9190, 'n_words': 2065, 'path': 'notebooks/04_analysis_1st-level.ipynb', 'repo_name': 'yarikoptic/fmriflows', 'size': '31852'}