--- a
+++ b/fmri_fsl_1stLvl.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python
+# %%
+"""
+Created on Wed Dec  4 14:29:06 2019
+
+@author: Or Duek
+1st level analysis using FSL output
+In this one we smooth using SUSAN, which takes longer. 
+"""
+# %%
+from __future__ import print_function
+from __future__ import division
+from builtins import str
+from builtins import range
+
+import os  # system functions
+
+import nipype.interfaces.io as nio  # Data i/o
+import nipype.interfaces.fsl as fsl  # fsl
+import nipype.interfaces.utility as util  # utility
+import nipype.pipeline.engine as pe  # pypeline engine
+import nipype.algorithms.modelgen as model  # model generation
+#import nipype.algorithms.rapidart as ra  # artifact detection
+#from nipype.workflows.fmri.fsl.preprocess import create_susan_smooth
+from nipype.interfaces.utility import Function
+"""
+The output file format for FSL routines is being set to compressed NIFTI.
+"""
+
+fsl.FSLCommand.set_default_output_type('NIFTI_GZ')
+
+
+data_dir = '/gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/CPT/BIDS/derivatives/'
+
+removeTR = 4
+fwhm = 4
+tr = 3
+
+output_dir = '/gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/'
+# %% Methods
+def _bids2nipypeinfo(in_file, events_file, regressors_file, removeTR = 4,
+                     regressors_names=None,
+                     motion_columns=None,
+                     decimals=3, amplitude=1.0):
+    from pathlib import Path
+    import numpy as np
+    import pandas as pd
+    from nipype.interfaces.base.support import Bunch
+    
+    # Process the events file
+    events = pd.read_csv(events_file, sep=',')
+
+    bunch_fields = ['onsets', 'durations', 'amplitudes']
+
+    if not motion_columns:
+        from itertools import product
+        motion_columns = ['_'.join(v) for v in product(('trans', 'rot'), 'xyz')]
+
+    out_motion = Path('motion.par').resolve()
+
+    regress_data = pd.read_csv(regressors_file, sep=r'\s+')
+    np.savetxt(out_motion, regress_data[motion_columns].values[removeTR:,], '%g')
+    if regressors_names is None:
+        regressors_names = sorted(set(regress_data.columns) - set(motion_columns))
+
+    if regressors_names:
+        bunch_fields += ['regressor_names']
+        bunch_fields += ['regressors']
+
+    runinfo = Bunch(
+        scans=in_file,
+        conditions=list(set(events.trial_type.values)),
+        **{k: [] for k in bunch_fields})
+
+    for condition in runinfo.conditions:
+        event = events[events.trial_type.str.match(condition)]
+
+        runinfo.onsets.append(np.round(event.onset.values-removeTR, 3).tolist()) # added -removeTR to align to the onsets after removing X number of TRs from the scan
+        runinfo.durations.append(np.round(event.duration.values, 3).tolist())
+        if 'amplitudes' in events.columns:
+            runinfo.amplitudes.append(np.round(event.amplitudes.values, 3).tolist())
+        else:
+            runinfo.amplitudes.append([amplitude] * len(event))
+
+    if 'regressor_names' in bunch_fields:
+        runinfo.regressor_names = regressors_names
+        runinfo.regressors = regress_data[regressors_names].fillna(0.0).values[removeTR:,].T.tolist() # adding removeTR to cut the first rows
+
+    return [runinfo], str(out_motion)
+# %%
+subject_list = ['2001', '2002','2004','2005','2008','2010','2012','2013','2015','2017','2021','2022','2022',
+'2023','2024','2025','2026','2027','2028','2032','2033','2034','2036','2037','2038','2039','2042',
+'2043','2044','2045','2047','2048','2050','2051','2052','2053','2054','2055','2056','2058','2059','2062',
+'2063','2064']
+# # Map field names to individual subject runs. 
+
+
+
+infosource = pe.Node(util.IdentityInterface(fields=['subject_id'
+                                            ],
+                                    ),
+                  name="infosource")
+infosource.iterables = [('subject_id', subject_list)]
+
+# SelectFiles - to grab the data (alternativ to DataGrabber)
+templates = {'func': data_dir +  '/fmriprep/sub-{subject_id}/func/sub-{subject_id}_task-imagery_space-MNI152NLin2009cAsym_res-2_desc-preproc_bold.nii.gz',
+             'mask': data_dir + '/fmriprep/sub-{subject_id}/func/sub-{subject_id}_task-imagery_space-MNI152NLin2009cAsym_res-2_desc-brain_mask.nii.gz',
+             'regressors': data_dir + '/fmriprep/sub-{subject_id}/func/sub-{subject_id}_task-imagery_desc-confounds_timeseries.tsv',
+             'events':  '/home/oad4/CPT_ML/CPT_event_file.csv'}
+selectfiles = pe.Node(nio.SelectFiles(templates,
+                               ),
+                   name="selectfiles")
+# %%
+
+# Extract motion parameters from regressors file
+runinfo = pe.Node(util.Function(
+    input_names=['in_file', 'events_file', 'regressors_file', 'regressors_names', 'removeTR'],
+    function=_bids2nipypeinfo, output_names=['info', 'realign_file']),
+    name='runinfo')
+runinfo.inputs.removeTR = removeTR
+
+# Set the column names to be used from the confounds file
+runinfo.inputs.regressors_names = ['dvars', 'framewise_displacement'] + \
+    ['a_comp_cor_%02d' % i for i in range(6)] + ['cosine%02d' % i for i in range(4)]
+# %%
+skip = pe.Node(interface=fsl.ExtractROI(), name = 'skip') 
+skip.inputs.t_min = removeTR
+skip.inputs.t_size = -1
+
+# %%
+susan =  pe.Node(interface=fsl.SUSAN(), name = 'susan') #create_susan_smooth()
+susan.inputs.fwhm = fwhm
+susan.inputs.brightness_threshold = 1000.0
+
+
+# %%
+modelfit = pe.Workflow(name='modelfit', base_dir= output_dir)
+
+
+modelspec = pe.Node(interface=model.SpecifyModel(),                  
+                    name="modelspec")
+
+modelspec.inputs.input_units = 'secs'
+modelspec.inputs.time_repetition = tr
+modelspec.inputs.high_pass_filter_cutoff= 120
+
+## Building contrasts
+level1design = pe.Node(interface=fsl.Level1Design(), name="level1design")
+cont1 = ['trauma>neutral', 'T', ['trauma1', 'neutral1'], [1, -1]]
+cont2 = ['trauma>think_trauma', 'T', ['trauma1', 'think_trauma1'], [1, -1]]
+cont3 = ['trauma>stop_trauma', 'T', ['trauma1', 'go_trauma1'], [1, -1]]
+
+contrasts = [cont1, cont2, cont3]
+
+
+level1design.inputs.interscan_interval = tr
+level1design.inputs.bases = {'dgamma': {'derivs': False}}
+level1design.inputs.contrasts = contrasts
+level1design.inputs.model_serial_correlations = True    
+
+
+modelgen = pe.Node(
+    interface=fsl.FEATModel(),
+    name='modelgen',
+    )
+mask =  pe.Node(interface= fsl.maths.ApplyMask(), name = 'mask')
+
+
+modelestimate = pe.Node(
+    interface=fsl.FILMGLS(smooth_autocorr=True, mask_size=5, threshold=1000),
+    name='modelestimate',
+    )
+
+
+# %%
+modelfit.connect([
+    (infosource, selectfiles, [('subject_id', 'subject_id')]),
+    (selectfiles, runinfo, [('events','events_file'),('regressors','regressors_file')]),
+    (selectfiles, skip,[('func','in_file')]),
+    (skip,susan,[('roi_file','in_file')]),
+    
+    (susan, runinfo, [('smoothed_file', 'in_file')]),
+    (susan, modelspec, [('smoothed_file', 'functional_runs')]),
+    (runinfo, modelspec, [('info', 'subject_info'), ('realign_file', 'realignment_parameters')]),
+    (modelspec, level1design, [('session_info', 'session_info')]),
+    (level1design, modelgen, [('fsf_files', 'fsf_file'), ('ev_files',
+                                                          'ev_files')]),
+   # (susan, changeTosrting, [('outputnode.smoothed_files', 'arr')]),
+    (susan, mask, [('smoothed_file', 'in_file')]),
+    (selectfiles, mask, [('mask', 'mask_file')]),
+    (mask, modelestimate, [('out_file','in_file')]),
+    (modelgen, modelestimate, [('design_file', 'design_file'),('con_file', 'tcon_file'),('fcon_file','fcon_file')]),
+    
+])
+# %%
+modelfit.run('MultiProc', plugin_args={'n_procs': 5})