=================================================================
= general outline of this complete analysis
=================================================================

# ======================================================================
# prepare input data
#    - in this case we have a BIDS tree format
#       - only hard part is creating reasonable timing files
#    - else
#       - for anat/EPI creation, consider: Dimon or dcm2niix_afni
#       - need mechanism to record and extract actual stimulus timing
#         at scanner
#       - it is VERY important to handle pre-steady state time
#          - adjust timing for number of volumes thrown out
#          - scanner might do this automatically, need to know
#          - never give stimulus during pre-SS time
#
# ==> now have anat, EPI, and stim timing files


# ------------------------------------------------------------
# no FreeSurfer in this case
#    - but if FS is used, need 1 mm^3 even-voxel grid


# ======================================================================
# skull-strip anat and compute non-linear warp to template
#    - run before any task analysis, apply to all (tasks, rest)
#    - see afni_scripts/p.1.sswarper
#          afni_scripts/c.ss.1.sswarper

   @SSwarper -input $indata_root/$sid/anat/${sid}_T1w.nii.gz  \
             -base MNI152_2009_template_SSW.nii.gz            \
             -subid $sid -odir .


# ======================================================================
# create timing files: ** this can be the hardest part **
#    - see afni_scripts/p.2.make.timing
#
# convert ${sid}_task-pamenc_events.tsv into timing
#    - TASK and CONTROL conditions
#    - using duration modulation, so attach duration to onset
#
#    4.52:0.65  19.51:0.72  49.52:0.63  70.52:1.85  ...


# ======================================================================
# run afni_proc.py
#    - see afni_scripts/p.3.run.AP
#          afni_scripts/c.ss.3.AP.pamenc

set task = pamenc
set SSdir = AFNI_01_SSwarper
set APdir = AFNI_02_${task}

# anat: input stripped anat, unifized anat, and WARP info
set anat_orig     = .../$sid/anatU.$sid.nii
set anat_ss       = .../$sid/anatSS.$sid.nii
set anat_std      = .../$sid/anatQQ.$sid.nii
set anat_warp_aff = .../$sid/anatQQ.$sid.aff12.1D
set anat_warp_nl  = .../$sid/anatQQ.${sid}_WARP.nii

# EPI: input just one run here
set EPI_files     = ( $BIDS/$sid/func/${sid}_task-${task}_bold.nii.gz )

# timing files were created above
set timing_files  = ( .../$APdir/$sid/timing/times.{CONTROL,TASK}.txt )
set stim_classes  = ( CONTROL TASK )


# ------------------------------
# task-pamenc_bold.json shows slice timing of alt+z2 (missing from nii.gz)
# blur in mask, and use higher 6 mm FWHM (voxels are 3x3x4)

afni_proc.py -subj_id $sid                                         \
   # data inputs and options specific to them 
      -copy_anat $anat_ss                                          \
      -anat_has_skull no                                           \
      -anat_follower anat_w_skull anat $anat_orig                  \

      -dsets $EPI_files                                            \
      -tcat_remove_first_trs 0                                     \
      -tshift_opts_ts -tpattern alt+z2                             \

      -regress_stim_times $timing_files                            \
      -regress_stim_labels $stim_classes                           \
      -regress_stim_types AM1                                      \
      -regress_basis_multi dmBLOCK                                 \

   # preprocessing options
      -blocks tshift align tlrc volreg mask blur scale regress     \
      -align_opts_aea -giant_move -cost lpc+ZZ (-check_flip)       \
      -tlrc_base MNI152_2009_template_SSW.nii.gz                   \
      -tlrc_NL_warp                                                \
      -tlrc_NL_warped_dsets $anat_std $anat_warp_aff $anat_warp_nl \
      -volreg_align_to MIN_OUTLIER                                 \
      -volreg_align_e2a                                            \
      -volreg_tlrc_warp                                            \
      -mask_epi_anat yes                                           \
      -blur_size 6                                                 \
      -blur_in_mask yes                                            \

    # regression options
      -regress_motion_per_run                                      \
      -regress_censor_motion 0.3                                   \
      -regress_censor_outliers 0.05                                \
      -regress_compute_fitts           (just saving RAM)           \
      -regress_fout no                 (keeping small for class)   \
      -regress_opts_3dD                                            \
         -jobs $njobs                                              \
      -regress_3dD_stop                                            \
      -regress_reml_exec                                           \
      -regress_est_blur_errts                                      \

    # non-automatic QC options
      -radial_correlate yes                                        \
      -align_opts_aea ... -check_flip                              \
      -regress_make_ideal_sum sum_ideal.1D                         \
      -regress_run_clustsim no                                     \
      -html_review_style pythonic


# ======================================================================
# post-analysis QC and group prep

# write spreadsheet of out.ss_review text files
gen_ss_review_table.py -tablefile QC/ss_review_table.xls \
        -infiles sub*/s*.results/out.ss*.txt

# check fields that should not vary across subjects
gen_ss_review_table.py -outlier_sep space                       \
        -report_outliers 'AFNI version' VARY                    \
        -report_outliers 'num regs of interest' VARY            \
        -report_outliers 'final voxel resolution' VARY          \
        -report_outliers 'num TRs per run' VARY                 \
        -infiles sub*/s*.results/out.ss*.txt                    \
        -write_outliers QC/outliers.a.VARY.txt

# determine subjects to drop
gen_ss_review_table.py -outlier_sep space                       \
        -report_outliers 'censor fraction' GE 0.15              \
        -report_outliers 'average censored motion' GE 0.1       \
        -report_outliers 'max censored displacement' GE 8       \
        -infiles sub*/s*.results/out.ss*.txt                    \
        -write_outliers QC/outliers.b.short.txt

set bad_subs = ( `awk '/sub-/ {print $1}' QC/outliers.b.short.txt` )

# ------------------------------
# get ACF (blur) means for Monte Carlo cluster simulations
grep -h ACF sub*/*.results/out.ss*.txt | awk -F: '{print $2}'   \
        | 3dTstat -mean -prefix - 1D:stdin\'

# ------------------------------
# generate 70% group mask and related
3dmask_tool -input sub*/*.results/mask_epi_anat*.HEAD \
            -prefix QC/group_mask.7 -frac 0.7
3dmask_tool -input sub*/*.results/mask_epi_anat*.HEAD \
            -prefix QC/group_mask.inter -frac 1.0

# mean of masks
3dMean -prefix QC/mask.mean sub*/*.results/mask_epi_anat*.HEAD

# ------------------------------
# "time" series of final anats and EPIs for registration comparision
3dTcat -prefix QC/all.EPI.vr.tcat     sub*/*.results/final_epi_vr*.HEAD
3dTcat -prefix QC/all.anat.final.tcat sub*/*.results/anat_final*.HEAD


# ======================================================================
# group analysis (see run.ggc.9.bipolar.T-C)
# (no mask initially)

# note subject lists: CONTROL, BIPOLAR, dropped (from either)
set part_file = participants.short.tsv
set subs_control = ( `awk '/CONTROL/ {print $1}' QC/$part_file` )
set subs_bipolar = ( `awk '/BIPOLAR/ {print $1}' QC/$part_file` )
set subs_drop    = ( `cat QC/outliers.c.drop.subs.txt` )

set drop_opt     = ( -dset_sid_omit_list $subs_drop )
set template     = MNI152_2009_template_SSW.nii.gz
set mask         = group_mask.7+tlrc

set beta1        = "TASK#0_Coef"
set beta2        = "CONTROL#0_Coef"
set setlab       = ( task control )
set label        = g_bipolar.paried.T-C
set gdir         = group_analysis.$tindex.$label

set tt_script    = run.tt.$tindex.$label

# list ALL subject datasets, then specify which to use/drop
gen_group_command.py -command 3dttest++                 \
        -write_script $tt_script                        \
        -dsets ../sub-*/*.results/stats.sub*REML+tlrc.HEAD \
        -dset_sid_list $subs_bipolar                    \
        $drop_opt                                       \
        -dsets ../sub-*/*.results/stats.sub*REML+tlrc.HEAD \
        -dset_sid_list $subs_bipolar                    \
        $drop_opt                                       \
        -subj_prefix sub-                               \
        -set_labels $setlab                             \
        -subs_betas "$beta1" "$beta2"                   \
        -verb 2                                         \
        -options -paired                                \
        |& tee out.ggc

tcsh -x $tt_script |& tee out.$tt_script


# ======================================================================
# clusterize

# ------------------------------------------------------------
# run cluster simulation (from QC directory)
mkdir files_ClustSim
set params = ( `cat out.ACF.means.1D` )
3dClustSim -both -mask group_mask.7+tlrc -acf $params \
        -cmd 3dClustSim.ACF.cmd -prefix files_ClustSim/ClustSim.ACF

# ------------------------------------------------------------
# check cluster size needed (for NN1 bi-sided clustering)
set clust_file = files_ClustSim/ClustSim.ACF.NN1_bisided.1D
1d_tool.py -infile $clust_file -csim_show_clustsize
set min_clust = `1d_tool.py -infile $clust_file -csim_show_clustsize -verb 0`

# ------------------------------------------------------------
# check p=0.001 t-stat threshold (symmetric, pos/neg)
# note that we actually have 25 subjects per group, so t-test DF = 24
# (check CONTROL and BIPOLAR results)
3dinfo group_analysis.1*/ttest*.HEAD |& grep statpar
3dinfo group_analysis.2*/ttest*.HEAD |& grep statpar

# check p=0.001 t-stat threshold (symmetric, pos/neg)
set t_thresh = `cdf -p2t fitt 0.001 24`

# ------------------------------------------------------------
# and finally clusterize

# table
3dClusterize -mask group_mask.7+tlrc -inset ttest++_result+tlrc \
   -idat 0 -ithr 1 -NN 1 -clust_nvox 17 -bisided -3.7454 3.7454

# mask statistical results
3dcalc -a ttest++_result+tlrc -b group_mask.7+tlrc -expr 'a*b'  \
       -prefix ttr_masked

# make ROI cluster mask
3dClusterize -mask group_mask.7+tlrc -inset ttest++_result+tlrc \
   -idat 0 -ithr 1 -NN 1 -clust_nvox 17 -bisided -3.7454 3.7454 \
   -pref_map ttr_clust_rois

# get cluster overlaps from dataset MNI_caez_ml_18+tlrc
# (first get name of atlas)
grep -B 1 MNI_caez_ml_18 ~/abin/AFNI_atlas_spaces.niml
whereami -omask ttr_clust_rois+tlrc -atlas CA_ML_18_MNI > CAEZ_overlaps.txt

3dROIstats -mask ttr_clust_rois+tlrc -nzvoxels ttest++_result+tlrc