#!/bin/tcsh -f
set VERSION = '$Id: mri_glmfit-sim,v 1.66 2016/11/16 20:52:22 greve Exp $';
setenv LANG en_US.UTF-8

set inputargs = ($argv);

set glmdir = ()
set nulltype = ();
set nsim = ()
set thresh = ()
set csdbase = ()
set simsign = ()
set cwpvalthresh = .05
set Overwrite = 0;
set DoSim = 1;
set LF = ();
set AllDoneFile = ();
set Seed = ();
set fwhmOverride = ();
set fwhmAdd = 0;
set UseUniformPDF = 0;
set UniformPDFMin = ();
set UniformPDFMax = ();
set DiagCluster = 0;
set PBNodeType = ();
set UseCache = 0;
set PermForce = 0;
set PermResid = 0; # use residual when running permutation, ter Braak method
set UseGRF = 0;
set volsubject = fsaverage;
set subjectOverride = ();
set Bonferroni = 0;
set ReportCentroid = 0;
set annot = aparc
set IsPVR = 0;
set DoClusterMean = 1
set OutputAnnot = 1
set DoGRFOCNAnat = 0;
set y = ()

set nJobs = 1;
set DoBackground = 0;
set DoPBSubmit = 0;
set DoPoll = 0;
set tSleepSec = 10;

set CacheDir = $FREESURFER_HOME/average/mult-comp-cor
set CacheLabel = cortex;

set tmpdir = ();
set CleanUp = 1;
set PrintHelp = 0;

if($#argv == 0) goto usage_exit;
set n = `echo $argv | grep -e --help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
endif
set n = `echo $argv | grep -e --version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

if(! -e $glmdir) then
  echo "ERROR: cannot find $glmdir"
  exit 1;
endif

set glmfitlog = $glmdir/mri_glmfit.log
if(! -e $glmfitlog) then
  echo "ERROR: cannot find $glmfitlog"
  exit 1;
endif

if($nulltype != perm) then
  set fwhmfile = $glmdir/fwhm.dat
  if(! -e $fwhmfile) then
    echo "ERROR: cannot find $fwhmfile"
    exit 1;
  endif
  # This may be overridden
  set fwhm = `cat $fwhmfile`;
else
  set fwhm = 0;
endif

set glmfitcwd = `cat $glmfitlog | awk '{if($1 == "cwd") print $2}'`
if(0 && ! -e $glmfitcwd) then
  echo "ERROR: cannot find $glmfitcwd"
  exit 1;
endif
set glmfitcwd = `pwd`;

set anattype = volume;
set subject = ();
set hemi = ();
set surf = "white";
set wls = ();

set glmfitcwd = `cat $glmfitlog | awk '{if($1 == "cwd") print $2}'`

# Go through the original mri_glmfit command-line
set glmfitcmd0 = `cat $glmfitlog | awk '{if($1 == "cmdline") print $0}'`
set glmfitcmd = ($glmfitcmd0);
echo $glmfitcmd
set gd2mtx = dods
set UseTable = 0;
set label = ();
while($#glmfitcmd)
  set flag = $glmfitcmd[1]; shift glmfitcmd;
  #echo $flag
  switch($flag)
  case "--surf"
  case "--surface"
    set subject = $glmfitcmd[1]; shift glmfitcmd;
    set hemi    = $glmfitcmd[1]; shift glmfitcmd;
    set anattype = surface;
    if($#glmfitcmd != 0) then
      set c = `echo $glmfitcmd[1] | cut -c 1-2`;
      if("$c" != "--") then
        set surf = $glmfitcmd[1]; shift glmfitcmd;
      endif
    endif
    breaksw
  case "--fsgd"
    shift glmfitcmd;
    if($#glmfitcmd > 0) then
      if($glmfitcmd[1] == doss) set gd2mtx = doss
    endif
    breaksw
  case "--label"
    set label = $glmfitcmd[1]; shift glmfitcmd;
    breaksw
  case "--perm-force"
    # Probably never gets here
    set PermForce = 1;
    breaksw
  case "--table"
  case "--y"
    if($#y == 0) set y = $glmfitcmd[1]; 
    shift glmfitcmd;
    if($DoClusterMean) then
      if(! -e $y) then
        set y = $glmfitcwd/$y
        if(! -e $y) then
          set y = $glmdir/`basename $y`
          if(! -e $y) then
            echo "ERROR: cannot find $y"
            echo "If you do not have this file, run with --no-y as the first option"
            exit 1;
          endif
        endif
      endif
    endif
    if("$flag" == "--table") set UseTable = 1;
    breaksw
  case "--wls"
    set wls = $glmfitcwd/$glmfitcmd[1]; shift glmfitcmd;
    if(! -e $wls) then
      set wls = $glmdir/`basename $wls`
      if(! -e $wls) then
        echo "ERROR: cannot find $wls"
        exit 1;
      endif
    endif
    breaksw
  # Ignore these options that have no args
  case "cmdline"
  case "mri_glmfit"
  case "--osgm"
  case "--prune"
  case "--no-prune"
  case "--pca"
  case "--synth"
  case "--allowsubjrep"
  case "--illcond"
  case "--debug"
  case "--synth"
  case "--cortex"
  case "--kurtosis"
  case "--nii"
  case "--nii.gz"
  case "--rescale-x"
  case "--no-rescale-x"
  case "--fisher"
  case "--no-pcc"
  case "dods"
  case "mri_glmfit.bin"
    breaksw
  # Ignore these options that have one arg
  case "--glmdir"
  case "--o"
  case "--C"
  case "--X"
  case "--w"
  case "--fwhm"
  case "--mask"
  case "--seed"
  case "--var-fwhm"
  case "--yffxvar"
  case "--ffxdof"
  case "--ffxdofdat"
  case "--exclude-frame"
    shift glmfitcmd;
    breaksw
  case "--pvr"
    set IsPVR = 1;
    shift glmfitcmd;
    breaksw
  # Ignore these options that have two args
  case "--rand-split"
    shift glmfitcmd;
    shift glmfitcmd;
    breaksw
  default
    echo ""
    echo "INFO: unrecognized mri_glmfit cmd option $flag"
    echo " ... probably not a problem"
    echo ""
    #sleep 10
    #exit 1;
    breaksw 
  endsw
end

set mask = ();
if(! $UseTable) then
  set mask = `stem2fname $glmdir/mask`
  if($status) then
    echo "$mask"
    exit 1;
  endif
endif


if(0 && $UseGRF && $anattype != volume) then
  # not true anymore
  echo "ERROR: can only use --grf with volume"
  exit 1;
endif
if($UseGRF) then
  set csdbase = "grf.th$thresh.$simsign"
  #set tmp = `echo "$thresh*10"|bc -l`
  #set tmp = `printf %02d $tmp`
  #set csdbase = "grf.th$tmp.$simsign"
endif

if($#subjectOverride) set subject = $subjectOverride

# Make sure the subject exists
if($#subject != 0) then
  if(! -e $SUBJECTS_DIR/$subject) then
    echo "ERROR: cannot find $subject in $SUBJECTS_DIR"
    exit 1;
  endif
  echo "SURFACE: $subject $hemi"
endif

# Get a list of contrasts, etc
set FSGDhasCon = 0;
if(-e $glmdir/y.fsgd) then
  # check whether contrasts were specified in the FSGD
  set FSGDhasCon = `cat $glmdir/y.fsgd | awk '{if($1 == "Contrast") print 1}' | wc -l`
endif
if($#tmpdir == 0) set tmpdir = $glmdir/tmp.mri_glmfit-sim-$$
mkdir -p $tmpdir
set clist = ($glmdir/*/C.dat)
if($status) then
  echo "ERROR: cannot find any contrasts"
  exit 1;
endif
set clist2 = ();
set conlist = ();
foreach c ($clist)
  set tmp = `dirname $c`
  set conname = `basename $tmp`
  set conlist = ($conlist $conname);
  set confile = $tmpdir/$conname.mtx
  cp $c $confile
  # If Cons in the FSGD, dont spec on glmfit cmd line
  if($FSGDhasCon == 0)  set clist2 = ($clist2 --C $confile);

  # Make sure sig is there
  set sig = `stem2fname $glmdir/$conname/sig`
  if($status) then
    echo "$sig"
    exit 1;
  endif

  # Make sure that csd is not there
  if($DoSim) then
    ls $glmdir/csd/$csdbase.j???-$conname.csd >& /dev/null
    set csdexists = (! $status)
    if($csdexists) then
      set csdfiles = ($glmdir/csd/$csdbase.j???-$conname.csd)
      if(! $Overwrite) then
        echo "ERROR: $csdfiles already exists"
        echo "Delete it or run with --overwrite"
        exit 1;
      endif
      rm -f $csdfiles
    endif
  endif
end

set StartDate = `date`

# Create log file
if($#LF == 0) then
  if($UseCache == 0) then
    set LF = $glmdir/$csdbase.mri_glmfit-sim.log
  else
    set LF = $glmdir/cache.mri_glmfit-sim.log
  endif
  rm -f $LF
endif
echo "log file is $LF"

echo "" | tee -a $LF
echo "cd `pwd`" | tee -a $LF
echo $0 | tee -a $LF
echo $inputargs | tee -a $LF
echo "" | tee -a $LF
echo $VERSION  | tee -a $LF
date | tee -a $LF
uname -a | tee -a $LF
echo $user | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "FREESURFER_HOME $FREESURFER_HOME" | tee -a $LF
echo "" | tee -a $LF
echo "Original mri_glmfit command line:" | tee -a $LF
echo $glmfitcmd0 | tee -a $LF
echo "" | tee -a $LF
echo "DoSim = $DoSim" | tee -a $LF
echo "UseCache = $UseCache" | tee -a $LF
echo "DoPoll = $DoPoll" | tee -a $LF
echo "DoPBSubmit = $DoPBSubmit" | tee -a $LF
echo "DoBackground = $DoBackground" | tee -a $LF
echo "DiagCluster = $DiagCluster" | tee -a $LF
echo "gd2mtx = $gd2mtx" | tee -a $LF
if($#Seed) echo "Seed = $Seed" | tee -a $LF
echo "fwhm = $fwhm" | tee -a $LF
if($#fwhmOverride) then 
  # Change fwhm from value in glmdir to override value
  echo "fwhmOverride = $fwhmOverride" | tee -a $LF
  set fwhm = $fwhmOverride
endif

if($UseCache) then
  #set fwhm = `cat $glmdir/fwhm.dat`;
  set fwhmStr = `perl -e "printf('"'%02d'"',int ( $fwhm+0.5+$fwhmAdd ) ) "`
  if($fwhmStr == 00)   set fwhmStr = 01;
  set threshStr = `perl -e "print 10*$thresh"`;
  set csdCache = $CacheDir/$subject/$hemi/$CacheLabel/fwhm$fwhmStr/$simsign/th$threshStr/mc-z.csd
  if(! -e $csdCache) then
    echo "ERROR: cannot find $csdCache" | tee -a $LF
    exit 1;
  endif
  echo CSD $csdCache | tee -a $LF
  if($#csdbase == 0) set csdbase = cache.th$threshStr.$simsign
endif

if($DoSim) then
  # Use fsgd if there, if not use X
  set fsgd = ();
  set X = ();
  if(-e $glmdir/y.fsgd) then
    set fsgd = $glmdir/y.fsgd
  else
    set X = $glmdir/Xg.dat
  endif

  # Number of sim iterations per job
  set nSimPerJob = `echo "$nsim/$nJobs" | bc`;
  echo "nSimPerJob = $nSimPerJob"| tee -a $LF

  if($DoPoll) then
    set PollDir =  $glmdir/csd/poll
    mkdir -p $PollDir
    # User can create this file to stop simulation
    set StopFile = $PollDir/StopSim
    rm -f $StopFile
  endif

  # ----------------- Run glmfit simulation ----------------------
  # This may take a while
  mkdir -p $glmdir/csd
  @ nthJob = 1
  while($nthJob <= $nJobs)
    set jpad = `printf j%03d $nthJob`
    echo "$nthJob/$nJobs `date`" | tee -a $LF

    set cmd = (mri_glmfit $clist2)
    set cmd = ($cmd --sim $nulltype $nSimPerJob $thresh \
       $glmdir/csd/$csdbase.$jpad)
    set cmd = ($cmd --sim-sign $simsign --fwhm $fwhm)
    if($#fsgd) set cmd = ($cmd --fsgd $fsgd $gd2mtx)
    if($#X)    set cmd = ($cmd --X $X)
    if($#wls)  set cmd = ($cmd --wls $wls)
    # --y must go after --fsgd
    if(! $UseTable) then
      set cmd = ($cmd --mask $mask)
      if($PermResid == 0) set cmd = ($cmd --y $y)
      if($PermResid == 1) set cmd = ($cmd --y $residual)
    endif
    if($UseTable)   set cmd = ($cmd --table $y)
    #if($#label) set cmd = ($cmd --label $label)
    if($PermForce) set cmd = ($cmd --perm-force)
    if($anattype == surface) set cmd = ($cmd --surf $subject $hemi $surf);
    if($#Seed) set cmd = ($cmd --seed $Seed);
    if($UseUniformPDF) set cmd = ($cmd --uniform $UniformPDFMin $UniformPDFMax);
    if($DiagCluster) set cmd = ($cmd --diag-cluster --debug)
    if($?verbose) set cmd = ($cmd --debug)
    if($DoPoll) then
      set DoneFile = $PollDir/done.$csdbase.$jpad
      rm -f $DoneFile
      set cmd = ($cmd --sim-done $DoneFile)
    endif

    echo $cmd | tee -a $LF
    if(! $DoPoll) then
      # Run directly in shell
      $cmd | tee -a $LF
      if($status && ! $DiagCluster) exit 1;
      date | tee -a $LF

      # ----------------------------------------------------
      if($DiagCluster) then
        set conname = $conlist[1];
        set sig = ./$conname-sig.mgh
        set csd = $glmdir/csd/$csdbase.$jpad-$conname.csd
        set sum = $glmdir/csd/$csdbase.$jpad-$conname.diag

        if($anattype == surface) then
          set cmd = (mri_surfcluster --subject $subject \
            --hemi $hemi --surf $surf --cwpvalthresh $cwpvalthresh)
          if($ReportCentroid) set cmd = ($cmd --centroid);
          set c = 4;
        endif
        if($anattype == volume) then
          set cmd = (mri_volcluster)
          set c = 3;
        endif
        set cmd = ($cmd --in $sig --mask $mask --cwpvalthresh $cwpvalthresh \
         --sum $sum --thmin $thresh --sign $simsign)
        echo $cmd | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;

        set maxclustB = `grep -v \# $csd | awk '{print $3}' | sort -nr | head -n 1`;
        set maxclustA = `grep -v \# $sum | awk -v c=$c '{print $c}' | sort -nr | head -n 1`;
        echo "MaxClusterSizeA $maxclustA"| tee -a $LF
        echo "MaxClusterSizeB $maxclustB"| tee -a $LF
        exit 1;
      endif

      # ----------------------------------------------------
    else
      if($DoBackground) then
        # Run in background
        $cmd >>& $LF &
      endif
      if($DoPBSubmit) then
        # pbsubmit
        echo "PBSUBMIT ---------------------------------" | tee -a $LF
        pbsubmit $PBNodeType -c "$cmd >> $LF" | tee -a $LF
        sleep 10;
        echo "" | tee -a $LF
      endif
    endif
    @ nthJob = $nthJob + 1
  end

  if($DoPoll) then
    # Possible bug: if job above dies without creating done file
    echo "Polling" | tee -a $LF
    @ nthPoll = 0;
    set Done = 0;
    while(! $Done)
      if(-e $StopFile) then
        # Something has gone wrong, maybe user deleted to halt
        echo "WARNING: simulation aborted by user" | tee -a $LF
        exit 1;
      endif
      if(! -e $PollDir) then
        # Something has gone wrong, maybe user deleted to halt
        echo "ERROR: cannot find $PollDir, aborting" | tee -a $LF
        exit 1;
      endif
      @ nthPoll = $nthPoll + 1;
      set Done = 1;
      @ nthJob = 1
      while($nthJob <= $nJobs)
        set jpad = `printf j%03d $nthJob`
        set DoneFile = $PollDir/done.$csdbase.$jpad
        if(! -e $DoneFile) then
          echo "Poll $nthPoll job $nthJob `date`" | tee -a $LF
          set Done = 0;
          break;
        endif
      @ nthJob = $nthJob + 1
      end
      sleep $tSleepSec 
    end
  endif

endif #DoSim

if($DoPBSubmit) then
  # pbsubmit the surfcluster part
  set cmd = (mri_glmfit-sim --glmdir $glmdir  --tmp $tmpdir)
  set cmd = ($cmd --sim-sign $simsign --no-sim $csdbase --cwpvalthresh $cwpvalthresh)
  set cmd = ($cmd --log $LF)
  echo "------------------------------------------------" | tee -a $LF
  echo $cmd  | tee -a $LF
  pbsubmit $PBNodeType -c "$cmd"
  sleep 5;
  echo "" 
  echo "Just submitted this as a job, and exiting now."
  echo ""
  # Could wait until it is done by polling for AllDoneFile
  exit 0;
endif

# -------------- now run clustering -------------------- #
foreach conname ($conlist)
  if($UseTable) continue;

  set csdlist = ();
  if(! $UseGRF) then
    if($UseCache) then
      set csdfiles = ($csdCache)
    else
      ls $glmdir/csd/$csdbase.j???-$conname.csd >& /dev/null
      if($status) then
        echo "ERROR: cannot find any csd files"| tee -a $LF
        exit 1; 
      endif
      set csdfiles = ($glmdir/csd/$csdbase.j???-$conname.csd)
    endif

    set csdlist = ();
    foreach csd ($csdfiles)
      set csdlist = ($csdlist --csd $csd);
    end
  endif

  set sig = `stem2fname $glmdir/$conname/sig`
  set ext = `fname2ext $sig`
  set vwsig = $glmdir/$conname/$csdbase.sig.voxel.$ext
  set vwsigmax = $glmdir/$conname/$csdbase.sig.voxel.max.dat
  set cwsig = $glmdir/$conname/$csdbase.sig.cluster.$ext
  set msig  = $glmdir/$conname/$csdbase.sig.masked.$ext
  set ocn   = $glmdir/$conname/$csdbase.sig.ocn.$ext
  set oannot = $glmdir/$conname/$csdbase.sig.ocn.annot
  set sum   = $glmdir/$conname/$csdbase.sig.cluster.summary
  set csdpdf = $glmdir/$conname/$csdbase.pdf.dat

  if($anattype == surface) then
    set cmd = (mri_surfcluster --in $sig  --mask $mask \
      --cwsig $cwsig --sum $sum --ocn $ocn \
      --annot $annot --cwpvalthresh $cwpvalthresh --o $msig --no-fixmni)
    if(! $UseGRF) set cmd = ($cmd $csdlist --csdpdf $csdpdf --vwsig $vwsig --vwsigmax $vwsigmax)
    if($UseGRF)   set cmd = ($cmd --fwhm $fwhm --hemi $hemi --subject $subject --thmin $thresh --sign $simsign)

    if($OutputAnnot) set cmd = ($cmd --oannot $oannot)
    if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
    if($ReportCentroid) set cmd = ($cmd --centroid);
    set cmd = ($cmd --surf $surf); # bug to have to do this
    echo $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
    if($DoClusterMean) then
      # Create an average "waveform" for each cluster
      set avgwf = $glmdir/$conname/$csdbase.y.ocn.dat
      set cmd = (mri_segstats --seg $ocn --exclude 0  \
       --i $y  --avgwf $avgwf --sum /tmp/mri_glmfit-sim.junk.$$)
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      rm -f /tmp/mri_glmfit-sim.junk.$$
    endif
  endif

  if($anattype == volume) then
    mri_info --o  /tmp/mri_glmfit-sim.junk.$$.cres --cres $sig
    set voxsize = `cat /tmp/mri_glmfit-sim.junk.$$.cres`
    if($voxsize != 2) then
      echo "ERROR: voxel size of input is $voxsize, expecting 2"
      exit 1
    endif
    set reg = $FREESURFER_HOME/subjects/$volsubject/mri.2mm/reg.2mm.dat 
    set aseg = aparc+aseg.mgz
    if(! -e $FREESURFER_HOME/subjects/$volsubject/mri/$aseg) set aseg = aseg.mgz
    set cmd = (mri_volcluster --in $sig --mask $mask --reg $reg --no-fixmni\
      --cwsig $cwsig --sum $sum --ocn $ocn --cwpvalthresh $cwpvalthresh \
      --seg $volsubject $aseg --out $msig --allowdiag)
    if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
    if(! $UseGRF) set cmd = ($cmd --csdpdf $csdpdf $csdlist --vwsig $vwsig  --vwsigmax $vwsigmax)
    if($UseGRF)   set cmd = ($cmd --fwhmdat $glmdir/fwhm.dat --sign $simsign --thmin $thresh  --sign $simsign)
    echo $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;

    if($DoClusterMean) then
      # Create an average "waveform" for each cluster
      set avgwf = $glmdir/$conname/$csdbase.y.ocn.dat
      set cmd = (mri_segstats --seg $ocn --exclude 0 \
       --i $y --avgwf $avgwf --sum /tmp/mri_glmfit-sim.junk.$$)
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      rm -f /tmp/mri_glmfit-sim.junk.$$
    endif

    if($UseGRF && $DoGRFOCNAnat) then
      # Map the OCN into anat space, takes a lot of disk space
      set ocnanat = $glmdir/$conname/$csdbase.sig.ocn.anat.$ext
      set cmd = (mri_vol2vol --interp nearest --mov $ocn --regheader \
        --targ $SUBJECTS_DIR/fsaverage/mri/orig.mgz \
        --o $ocnanat --no-save-reg)
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
    endif

  endif

  # For some reason cwsig has 4 frames, just take 1st
  set cmd = (mri_convert $cwsig $cwsig --frame 0)
  echo $cmd | tee -a $LF
  $cmd  | tee -a $LF
  if($status) exit 1;
  
end

if($CleanUp) rm -r $tmpdir

if($#AllDoneFile) touch $AllDoneFile

echo $StartDate | tee -a $LF
date | tee -a $LF
echo "mri_glmfit-sim done" | tee -a $LF

exit 0
###############################################

############--------------##################
parse_args:
set cmdline = ($argv);
while( $#argv != 0 )

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "--glmdir":
      if( $#argv < 0) goto arg1err;
      set glmdir = $argv[1]; shift;
      breaksw

    case "--y":
      if($#argv < 0) goto arg1err;
      # Override input in fsgd file
      set y = $argv[1]; shift;
      if(! -e $y) then
        echo "ERROR: cannot find $y"
        exit 1
      endif
      breaksw

    case "--mczsim":
    case "--mcsim":
    case "--cache":
      if( $#argv < 2) goto arg2err;
      set thresh = $argv[1]; shift;
      set simsign = $argv[1]; shift;
      set UseCache = 1;
      set DoSim = 0;
      set thresh = `printf %2.1f $thresh`
      if($thresh != 1.3 && $thresh != 2.0 && $thresh != 2.3 && \
         $thresh != 3.0 && $thresh != 3.3 && $thresh != 4.0 && $thresh != 1.6) then
         echo "ERROR: thresh = $thresh, must be 1.3, 2.0, 2.3, 3.0, 3.3, 4.0, 1.6"
         exit 1;
      endif
      breaksw

    case "--mczsim-dir":
    case "--cache-dir":
      if( $#argv < 1) goto arg1err;
      set CacheDir = $argv[1]; shift;
      set UseCache = 1;
      set DoSim = 0;
      breaksw

    case "--mczsim-label":
    case "--cache-label":
      if( $#argv < 1) goto arg1err;
      set CacheLabel = $argv[1]; shift;
      set UseCache = 1;
      set DoSim = 0;
      breaksw

    case "--vol-subject":
      if( $#argv < 1) goto arg1err;
      set volsubject = $argv[1]; shift;
      breaksw

    case "--grf":
      if( $#argv < 2) goto arg2err;
      set thresh = $argv[1]; shift;
      set simsign = $argv[1]; shift;
      #if($simsign != pos && $simsign != neg) then
        #echo "ERROR: must use pos or neg with --grf"
        #exit 1;
      #endif
      set UseGRF = 1;
      set DoSim = 0;
      breaksw

    case "--perm":
      if( $#argv < 3) goto arg4err;
      set nulltype = perm;
      set nsim = $argv[1]; shift;
      set thresh = $argv[1]; shift;
      set simsign = $argv[1]; shift;
      set threshStr = `perl -e "print 10*$thresh"`;
      set csdbase = perm.th$threshStr.$simsign
      breaksw

    case "--sim":
      if( $#argv < 4) goto arg4err;
      set nulltype = $argv[1]; shift;
      set nsim = $argv[1]; shift;
      set thresh = $argv[1]; shift;
      set csdbase = $argv[1]; shift;
      breaksw

    case "--sim-sign":
      if( $#argv < 1) goto arg1err;
      set simsign = $argv[1]; shift;
      breaksw

    case "--subject-override":
      if( $#argv < 1) goto arg1err;
      set subjectOverride = $argv[1]; shift;
      breaksw

    case "--centroid":
      set ReportCentroid = 1;
      breaksw

    case "--log":
      if( $#argv < 1) goto arg1err;
      set LF = $argv[1]; shift;
      breaksw

    case "--a2009s":
      set annot = "aparc.a2009s"
      breaksw
    case "--annot":
      if( $#argv < 1) goto arg1err;
      set annot = $argv[1]; shift;
      breaksw

    case "--overwrite":
      set Overwrite = 1;
      breaksw

    case "--pbsubmit":
      if($#argv < 1) goto arg1err;
      set nJobs = $argv[1]; shift;
      set DoPBSubmit = 1;
      set DoPoll = 1;
      set DoBackground = 0;
      if(`hostname` != seychelles && `hostname` != launchpad) then
        hostname
        echo "ERROR: must be on seychelles or launchpad to pbsubmit"
        exit 1;
      endif
      if(`hostname` == seychelles) set PBNodeType = "-l nodes=1:opteron";
      breaksw

    case "--bg":
      if($#argv < 1) goto arg1err;
      set nJobs = $argv[1]; shift;
      set DoBackground = 1;
      set DoPoll = 1;
      set DoPBSubmit = 0;
      breaksw

    case "--no-bg":
      set DoBackground = 0;
      set DoPoll = 0;
      breaksw

    case "--sleep":
      if($#argv < 1) goto arg1err;
      set tSleepSec = $argv[1]; shift;
      breaksw

    case "--seed":
      if($#argv < 1) goto arg1err;
      set Seed = $argv[1]; shift;
      breaksw

    case "--fwhm-override":
      if($#argv < 1) goto arg1err;
      set fwhmOverride = $argv[1]; shift;
      breaksw

    case "--fwhm-add":
      if($#argv < 1) goto arg1err;
      set fwhmAdd = $argv[1]; shift;
      breaksw

    case "--perm-force":
      set PermForce = 1;
      breaksw
    case "--perm-resid":
      set PermResid = 1;
      set PermForce = 1; # purpose of using residual is to allow non-orthog
      breaksw

    case "--tmp":
      if($#argv < 1) goto arg1err;
      set tmpdir = $argv[1]; shift;
      breaksw

    case "--bonferroni":
      if($#argv < 1) goto arg1err;
      set Bonferroni = $argv[1]; shift;
      breaksw
    case "--2spaces":
      set Bonferroni = 2;
      breaksw
    case "--3spaces":
      set Bonferroni = 3;
      breaksw

    case "--cwpvalthresh":
    case "--cwp":
      if($#argv < 1) goto arg1err;
      set cwpvalthresh = $argv[1]; shift;
      breaksw

    case "--uniform":
      if($#argv < 2) goto arg2err;
      set UniformPDFMin = $argv[1]; shift;
      set UniformPDFMax = $argv[1]; shift;
      set UseUniformPDF = 1;
      breaksw

    case "--base":
      if($#argv < 1) goto arg1err;
      set csdbase = $argv[1]; shift;
      breaksw

    case "--no-sim":
      if($#argv < 1) goto arg1err;
      set csdbase = $argv[1]; shift;
      set DoSim = 0;
      breaksw

    case "--diag-cluster":
      set DiagCluster = 1;
      breaksw

    case "--no-out-annot":
      set OutputAnnot = 0
      breaksw

    case "--no-cluster-mean":
    case "--no-y":
      set DoClusterMean = 0;
      breaksw

    case "--no-grf-ocn-anat":
      set DoGRFOCNAnat = 0;
      breaksw
    case "--grf-ocn-anat":
      set DoGRFOCNAnat = 1;
      breaksw

    case "--debug":
      set verbose = 1;
      set echo = 1;
      breaksw

    default:
      echo ERROR: Flag $flag unrecognized. 
      echo $cmdline
      exit 1
      breaksw
  endsw

end

goto parse_args_return;
############--------------##################

############--------------##################
check_params:
if($#glmdir == 0) then
  echo "ERROR: must spec --glmdir"
  exit 1;
endif
if($DoSim) then
  if($#nulltype == 0) then
    echo "ERROR: must spec --sim"
    exit 1;
  endif
  if($#simsign == 0) then
    echo "ERROR: must spec --sim-sign"
    exit 1;
  endif
endif
set residual = ()
if($PermResid) then
  set stem = $glmdir/eres
  set residual = `stem2fname $stem`
  if($status) then
    echo "ERROR: cannot find residual $stem needed with --perm-resid"
    echo "  when running mri_glmfit, make sure to include --eres-save"
    exit 1;
  endif
endif
if($DiagCluster && $DoBackground) then
  echo "ERROR: cannot background with --diag-cluster"
  exit 1;
endif
if($DiagCluster && $DoPBSubmit) then
  echo "ERROR: cannot pbsubmit with --diag-cluster"
  exit 1;
endif
if($DiagCluster && $DoSim == 0) then
  echo "ERROR: must simulate with --diag-cluster"
  exit 1;
endif

goto check_params_return;
############--------------##################

############--------------##################
arg1err:
  echo "ERROR: flag $flag requires one argument"
  exit 1
############--------------##################

############--------------##################
arg2err:
  echo "ERROR: flag $flag requires two arguments"
  exit 1
############--------------##################

############--------------##################
arg4err:
  echo "ERROR: flag $flag requires four arguments"
  exit 1
############--------------##################

############--------------##################
usage_exit:
  echo ""
  echo "mri_glmfit-sim --help"
  echo ""
  echo " --glmdir glmdir"
  echo " --cwp cwthresh : cluster-wise p-value threshold ($cwpvalthresh)"
  echo ""
  echo "Use pre-computed z-based monte-carlo simulations (see mri_mcsim)"
  echo " --mczsim vwthreshold sign"
  echo "If you have run your own simulations (mri_mcsim), then specify "
  echo "   --mczsim-dir dir : default is FREESURFER_HOME/average/mult-comp-cor"
  echo "   --mczsim-label label : default is cortex"
  echo ""
  echo "Use permutation simulation"
  echo " --perm nsim vwthreshold sign"
  echo "   --perm-resid : add this with non-orthogonal designs to run permutation"
  echo "      on the residual rather than on the raw data (must have run mri_glmfit"
  echo "      with --eres-save)"
  echo ""
  echo "Use Gaussian Random Fields (GRF), volumes only"
  echo " --grf vwthreshold sign : -log10(p) and sign (pos, neg, abs)"
  echo ""
  echo "Correct for multiple spaces"
  echo " --2spaces : additional Bonferroni correction across 2 spaces (eg, lh, rh)"
  echo " --3spaces : additional Bonferroni correction across 3 spaces (eg, lh, rh, mni305)"
  echo ""
  echo " --overwrite : delete previous CSDs"
  echo " --bg njobs : divide sim into njobs and put in background"
  echo " --sleep tSleepSec : number of seconds to sleep between background polls "
  echo ""
  echo " --a2009s : use aparc.a2009s instead of aparc for reporting region of vertex max"
  echo " --annot annot : use annot instead of aparc for reporting region of vertex max"
  echo ""
  echo " --log logfile : default is csdbase.mri_glmfit-sim.log"
  echo ""
  echo " --base csdbase : override csdbase name"
  echo " --no-sim csdbase : do not simulate, only run cluster"
  echo " --seed seed : set simluation seed"
  echo " --fwhm-override fwhm : override fwhm in glmdir":
  echo " --fwhm-add fwhmAdd : add fwhmAdd to the estimated fwhm before addressing table":
  echo " --uniform min max : use uniform PDF instead of gaussian":
  echo " --no-out-annot : do not output a cluster annotation"
  echo " --no-cluster-mean : do not compute means of each subject in each cluster (was --no-y)"
  echo " --y yFile : use yFile instead of determining yFile from glmdir"
  echo " --centroid : report the coordinates/annotation of the centroid instead of max"
  echo " --help"
  echo ""

  if(! $PrintHelp) exit 1;

  echo $VERSION

  cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'

exit 1;

#---- Everything below here is printed out as part of help -----#
BEGINHELP

This program is a script to run corrections for multiple comparisons
(CMC) on the volume or the surface. It is a front-end for the
mri_glmfit, mri_surfcluster, and mri_volcluster programs.

To use this program:

1. Run mri_glmfit on your data, giving it the glm directory
   (--glmdir), and all the other arguments.

2. Choose a voxel-wise (cluster-forming) threshold vwthresh. The
   vwthresh takes the value -log10(p) where p is the p-value.  So if
   vwthresh = 2, then p<.01. You should understand the difference
   between the cluster-forming threshold and the cluster-wise
   threshold (--cwp). For mczsim, vwthresh must be 1.3, 2.0, 2.3, 3.0,
   3.3, or 4.0 (unless you have simulated more values). There are no
   good methods to choose this value. Typically 2 or 3 is used.
   It is not recommended to use 1.3 unless performing a permutation
   operation. 

3. Choose the sign of the effect, either positive (pos), negative
   (neg) or absolute/unsigned (abs). If you have an a priori
   hypothesis (eg, group1 > group2), then choose positive. If you do
   not expect a particular sign, then use abs. Choosing a sign will
   yield stronger results (ie, smaller p-values), but you will lose
   voxels with the opposite sign. NOTE: the "abs" for GRF is only
   approximate. For volumes, it may give very conservative results.

4. Choose the method. This will be either: (1) --mczsim z-based monte
   carlo simulation (See surfer.nmr.mgh.harvard.edu/fswiki/BuildYourOwnMonteCarlo).
   (2) --perm for permutation simulation, or (3) --grf for Gaussian Random Field
   method (analytical). For permution, you must choose the number of simulation
   iterations (typically 5000-1000); permutation may take a few hours to 
   complete. GRF is for volumes only.

Now run the command, eg,

mri_glmfit-sim --glmdir glmdir --mczsim 2 abs

or

mri_glmfit-sim --glmdir glmdir --sim perm 1000 3 abs

The output will be in each contrast folder and will be something like
the following (where csdbase will be "cache" with --mczsim or "perm"
with --perm or "grf" with --grf)

csdbase.sig.voxel.mgh 
csdbase.sig.cluster.mgh 
csdbase.sig.cluster.summary
csdbase.y.ocn.dat
csdbase.sig.ocn.mgh
csdbase.sig.ocn.annot (for surfaces only)
csdbase.sig.masked.mgh 

csdbase.sig.voxel.mgh - the sig volume corrected for multiple
comparisons on a voxel-wise basis. The threshhold and sign are
irrelevant. The value at each voxel is the corrected -log10(p-value)
for that voxel.

csdbase.sig.cluster.mgh - the sig volume corrected for multiple
comparisons on a cluster-wise basis. The value at each voxel
is the -log10(p), where p is the pvalue of the cluster at 
that voxel. If that voxel does not belong to a cluster, its
value will be 0. 

csdbase.sig.cluster.summary - this is the cluster summary table 
(a simple text file).

csdbase.y.ocn.dat - this is a summary of the input (y) over each
cluster. It has a column for each cluster. Each row is a subject. The
value is the average of the input (y) in each cluster. This is a
simple text file.

csdbase.sig.ocn.mgh - output cluster number volume. The value of the
voxel is the integer cluster number that that voxel belongs to. This
can be used like a segmentation (eg, to load into tkmedit or
use for running mri_segstats).

csdbase.sig.ocn.annot (for surfaces only) - this is an annotation
with the annotation name being the cluster number. This can be
used like any other annotation (eg, to load into tksurfer or 
as input for mri_segstats or mris_anatomical_stats).

csdbase.sig.masked.mgh - the original sig volume masked to show
only clusters.

OTHER OPTIONS

--overwrite

Overwrite existing CSD files when performing simulation. Default 
is to not overwrite so that you do not clobber CSD files that you
spent a weekend creating.

--bg njobs

CAREFUL!!!! WARNING!!! Divide the number of simulation interations
into njobs background jobs. If the simulation jobs die uncleanly for
some reason, this will never quit because it looks for a file created
by the simulation job before continuing. If this file is not created,
then it will keep polling for this file forever. ALSO: if you are not
running on a computer with multiple nodes, then this will not make it
run any faster (in fact, it will run slower). You can stop
backgrounded jobs creating a file called glmdir/csd/poll/StopSim. If
you are in the Martinos center, do not attempt to use --bg on the
launchpad computational cluster.

--cwp cwpvalthresh
--cwpvalthresh cwpvalthresh

Only report clusters that have p < cwpvalthresh. Default is .05. If
you want to change this after running a permutation simulation, you do
not need to re-run the entire simulation. Just run mri_glmfit with the
--no-sim option.

--2spaces
--3spaces

This adds an addition correction for multiple comparisons for when the
current correction is part of a larger analysis. For example, if you
are doing a structural study looking at both the left and the right
hemispheres, mri_glmfit-sim must be run separately for each
hemisphere. By default, it will give you p-values corrected only for
one hemisphere. To correct for both hemispheres use --2spaces. For an
fMRI study in which you are using both hemispheres and subcortical
structures, then use --3spaces. The correction applied is Bonferroni
(N=2 for --2spaces and N=3 for --3spaces).

--no-sim csdbase

Run this if you have already run the permutation stimulation but just
want to run the correction for multiple comparisons again to change
the clusterwise threshold. "csdbase" should be "perm"

--seed seed

Set the seed for simulation random number generator. This is mostly
used for debugging.

--fwhm-override fwhm 

Override fwhm in glmdir with given value. This is mostly
used for debugging.

--no-cluster-mean

Do not compute means of each subject in each cluster. This computation
requires the y input file which might not be available in all cases.
When using this option, put it first on the command-line.

--y yFile

Specify the GLM input (ie, argument to --y in mri_glmfit) rather
than determining it from the glmdir. This can be handy if the
y file was deleted at some point.

--perm-resid

This is used for non-orthogonal designs (eg, when there is a covariate
present). It runs the permutation using the residual of the analysis
as the input rather than the raw data (the residual will have the
covariate regressed out). You must have originally run mri_glmfit with
the --eres-save option. This method is an APPROXIMATION known as the
ter Braak method (ter Braak, 1992) that allows for permutation to be
run on non-orthogonal designs. See also Winkler, 2014.

ter Braak, C.J.F., 1992. Permutation versus bootstrap significance
tests in multiple regression and ANOVA. In: Jöckel, K.-H., Rothe, G.,
Sendler, W. (Eds.), Bootstrapping and Related
Techniques. No. 1989. Springer-Verlag, Berlin, pp. 79–86.

Winkler, et al, 2014. Permutation inference for the general linear
model. Neuroimage 92 (2014) 381–397.