#!/bin/tcsh -f
set VERSION = '$Id: mri_glmfit-sim,v 1.36.2.6 2013/06/05 17:06:31 greve Exp $';

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 UseUniformPDF = 0;
set UniformPDFMin = ();
set UniformPDFMax = ();
set DiagCluster = 0;
set PBNodeType = ();
set UseCache = 0;
set PermForce = 0;
set UseGRF = 0;
set volsubject = fsaverage;
set subjectOverride = ();
set Bonferroni = 0;
set ReportCentroid = 0;
set annot = aparc

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

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

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

if($nulltype != perm) then
  set fwhmfile = $glmdir/fwhm.dat
  if(! -e $fwhmfile) then
    echo "ERROR: cannot find $fwhm"
    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(! -e $glmfitcwd) then
  echo "ERROR: cannot find $glmfitcwd"
  exit 1;
endif

set anattype = volume;
set subject = ();
set hemi = ();
set surf = "white";
set y = ();
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 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 "--y"
    set y = $glmfitcwd/$glmfitcmd[1]; shift glmfitcmd;
    if(! -e $y) then
      set y = $glmdir/`basename $y`
      if(! -e $y) then
        echo "ERROR: cannot find $y"
        exit 1;
      endif
    endif
    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 "dods"
    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"
    shift glmfitcmd;
    breaksw
  default
    echo ""
    echo "WARNING: unrecognized mri_glmfit cmd option $flag"
    echo ""
    sleep 10
    #exit 1;
    breaksw 
  endsw
end

if($UseGRF && $anattype != volume) then
  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
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
  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+.5 ) ) "`
  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
  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 --y $y $clist2 --mask $mask)
    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)
    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)
    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)

  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 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 $csdlist --mask $mask \
      --cwsig $cwsig --vwsig $vwsig --sum $sum --ocn $ocn \
      --oannot $oannot --annot $annot --csdpdf $csdpdf \
      --cwpvalthresh $cwpvalthresh --o $msig --no-fixmni)
    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;
    # 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($anattype == volume) then
    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)
    if($Bonferroni) set cmd = ($cmd --bonferroni $Bonferroni)
    if(! $UseGRF) set cmd = ($cmd --csdpdf $csdpdf $csdlist --vwsig $vwsig )
    if($UseGRF)   set cmd = ($cmd --fwhmdat $glmdir/fwhm.dat --sign $simsign --thmin $thresh)
    echo $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;

    # 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.$$

    if($UseGRF) then
      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
date
echo "mri_glmfit-sim done"

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 "--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) then
         echo "ERROR: thresh = $thresh, must be 1.3, 2.0, 2.3, 3.0, 3.3, 4.0"
         exit 1;
      endif
      breaksw

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

    case "--cache-label":
      if( $#argv < 1) goto arg1err;
      set CacheLabel = $argv[1]; shift;
      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 "--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 "--perm-force":
      set PermForce = 1;
      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 "--no-sim":
      if($#argv < 1) goto arg1err;
      set csdbase = $argv[1]; shift;
      set DoSim = 0;
      breaksw

    case "--diag-cluster":
      set DiagCluster = 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

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"
  echo ""
  echo " --glmdir glmdir"
  echo " --cwp thresh : cluster p-value thresh ($cwpvalthresh)"
  echo ""
  echo " --cache threshold sign"
  echo "   --cache-dir dir : default is FREESURFER_HOME/average/mult-comp-cor"
  echo "   --cache-label label : default is cortex"
  echo ""
  echo " --grf threshold sign : -log10(p) and pos or neg (volumes only)"
  echo ""
  echo " --sim nulltype nsim threshold csdbase"
  echo " --sim-sign sign : abs,pos,neg"
  echo ""
  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 " --no-sim csdbase : do not simulate, only run cluster"
  echo " --seed seed : set simluation seed"
  echo " --fwhm-override fwhm : override fwhm in glmdir":
  echo " --uniform min max : use uniform PDF instead of gaussian":
  echo ""
  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 GLM-based simulations for corrections
for multiple comparisons 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 whether you want to use a pre-cached simulation or whether
   you want to run a new simulation. You should use a pre-cached
   simulation if: this is a surface-based analysis, you are using
   fsaverage (or equivalent), you have whole-brain data, you want to
   use a voxel-wise threshold of 1.3, 2.0, 2.3, 3.0, 3.3, or 4.0, and
   you are going to use an mc-z simulation.  Pre-cached will finish
   after a few seconds, your own simulation will take hours or days.

2.a If you want to use a pre-cached, then specify "--cache threshold sign"
    where threshold is the vertex-wise threshhold and sign is abs, pos,
    or neg. Note: the csdbase will be set to cache.thresh.simsign (eg,
    if the threshold was set to 1.3 and sign to pos, then cache.th13.pos.
    Skip the remaining steps below.

2.b If you want to run your own simulation, then choose the type of
    simulation (nulltype) you want to run. The choices are:
    mc-z - synthesized, smoothed z-field (fast)
    mc-full - fully synthesized, smoothed gaussian field (slowest)
    perm - permutation on original data (probably fasted one)

3. Choose the voxel/vertex-wise threshhold used to define clusters. This
   will be -log10(p), where p is the p-value. Eg, to use all voxels/vertices
   with p < .01, then set the threshhold = 2.

4. Choose the number of iterations (nsim). This is usually 10,000
   (which can take quite a while).

5. Select the CSD base name (csdbase). CSD stands for "Cluster
   Simulation Data". There will be a text file for each contrast with
   the information about the null distribution. This file will be
   called csdbase-contrast.csd, so the csdbase should be something
   descriptive. Eg, if you have mc-z, with threshold=2, and sign=pos,
   then setting csdbase to mc-z.pos.2 would be appropriate. If the 
   CSD file exists, then you must delete it or run with --overwrite.

6. Select the threshold sign (--sim-sign). This is the sign of the
   threshold used to create the clusters for cluster-wise correction
   of multiple comparisons. Options are abs (unsigned), pos
   (positive), and neg (negative). If you have an a priori hypothesis
   (eg, group1 > group2), then choose positive. If you don not know
   the sign, then use abs. Choosing a sign will yield stronger results
   (ie, smaller p-values), but you will lose the other sign.

Now run the command:

mri_glmfit-sim --glmdir glmdir \
  --cache threshold sign

or

mri_glmfit-sim --glmdir glmdir \
  --sim nulltype nsim threshold csdbase \
  --sim-sign sign

Example,

mri_glmfit-sim --glmdir lh.thickness.sm05 \
  --cache 3  --sim-sign pos

or

mri_glmfit-sim --glmdir lh.thickness.sm05 \
  --sim mc-z 10000 3 mc-z.pos.3 --sim-sign pos

This will run mri_glmfit in simulation mode. It will find all the
contrasts that you ran in the original invocation of mri_glmfit and
create CSD files for them. It will use the mask and, for mc
simulations, the FWHM found in the glmdir.

It will then run mri_surfcluster or mri_volcluster for
each contrast, using the sig volume as input and creating the
following files in each contrast directory:

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 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 stimulation but just want to run
the correction for multiple comparisons again to change the
clusterwise threshold. "csdbase" is the same that was used in the
orignial invocation to mri_glmfit-sim.  Eg, if this was your original
mri_glmfit-sim command:

  mri_glmfit-sim --glmdir lh.thickness.sm05 --cwp .05\
    --sim mc-z 10000 3 mc-z.pos.3 --sim-sign pos

But then then you wanted to re-run with a higher CWP, then

  mri_glmfit-sim --glmdir lh.thickness.sm05 --no-sim mc-z.pos.3 --cwp .1

--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.