#!/bin/tcsh -f
# groupstats

set VERSION = 'groupstats 8.2.0';

set fsgd = ();
set slist = ();
set hemilist = ();
set DoMaps = 1;
set tmpdir = ();
set cleanup = 1;
set LF = ();
set outdir = ()
set fwhmlist = ();
set surflist = (white)
#set aparclist = (aparc BA aparc.a2009s aparc.DKTatlas40)
set aparclist = (aparc BA aparc.a2009s)
set ASegLUT = $FREESURFER_HOME/ASegStatsLUT.txt
set WMParcLUT = $FREESURFER_HOME/WMParcStatsLUT.txt
set OverWrite = 0;
set mapmeaslist = (thickness area volume curv sulc w-g.pct) 
set DoAparcStats = 1
set DoAsegStats = 1
set DoWMParcStats = 1
set SrcSurfReg = sphere.reg
set Replace53 = "--replace53"
set Wait = 0
set tpollsec = -1

set inputargs = ($argv);
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:

set StartTime = `date`;
set tSecStart = `date '+%s'`;

mkdir -p $outdir/rois $outdir/maps
pushd $outdir > /dev/null
set outdir = `pwd`;
popd > /dev/null

# Set up log file
if($#LF == 0) set LF = $outdir/groupstats.log
if($LF != /dev/null) rm -f $LF
echo "Log file for groupstats" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
echo "" | tee -a $LF
echo "FREESURFER_HOME $FREESURFER_HOME" | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
uname -a  | tee -a $LF

if($#fsgd) then
  cp $fsgd $outdir/group.fsgd
else
  # Create FSGD file for OSGM
  set fsgd = $outdir/group.fsgd
  rm -f $fsgd
  echo "GroupDescriptorFile 1" >> $fsgd
  echo "Class osgm" >> $fsgd
  echo "Contrast osgm 1" >> $fsgd
  foreach s ($slist)
    echo "Input $s osgm" >> $fsgd
  end
endif

set bfile = $outdir/build-stamp.txt
rm -f $bfile
set b1 = ()
foreach s ($slist)
  set b = `cat $SUBJECTS_DIR/$s/scripts/build-stamp.txt`
  echo "$s $b" >> $bfile
  if($#b1 == 0) set b1 = $b
  if($b != $b1) then
    echo ""
    echo "WARNING: multiple builds detected, see $bfile" | tee -a $LF
    echo ""
    sleep 5
  endif
end

echo $mapmeaslist > $outdir/mapmeas.list.txt
echo $hemilist    > $outdir/hemi.list.txt
echo $aparclist   > $outdir/aparc.list.txt
echo $fwhmlist    > $outdir/fwhm.list.txt
echo $slist       > $outdir/subjectslist.txt

#========================================================

pushd $SUBJECTS_DIR
if(! -e fsaverage) then
  ln -s $FREESURFER/subjects/fsaverage
endif
popd

# BA label files changed name with version 6 to add "_exvivo"
# Detect it with this. Below, make output files use original name
set BAex = ""
set f0 = $SUBJECTS_DIR/$slist[1]/stats/lh.BA.stats
if(! -e $f0) then
  set f = $SUBJECTS_DIR/$slist[1]/stats/lh.BA_exvivo.stats
  if(! -e $f) then
    echo "ERROR: cannot find $f0 or $f"
    exit 1;
  endif
  set BAex = "_exvivo"
endif

# APARC Stats ----------------------------------------------------
if($DoAparcStats) then
foreach aparc ($aparclist)
  foreach surf ($surflist)
    if($aparc == aparc.a2009s && $surf == pial) continue;
    foreach hemi ($hemilist)
      foreach meas ($mapmeaslist)
        if($meas == curv || $meas == sulc) continue # curv breaks aparcstats2table
        if($meas == w-g.pct && ($aparc != aparc || $surf != white)) continue 
        set outfile = $outdir/rois/$aparc.$hemi.$surf.$meas.dat
        if($meas != w-g.pct) then
          set parc = $aparc
          if($aparc == BA) set parc = BA$BAex
          if($surf == pial) set parc = pial.$aparc; # only aparc
          if(! -e $outfile || $OverWrite) then
            set cmd = (aparcstats2table --subjects $slist --hemi $hemi \
               -m $meas -t $outfile -p $parc)
            echo "#C# $cmd" | tee -a $LF
            $cmd | tee -a $LF
            if($status) exit 1;
            set glmdir = $outdir/rois/glm.$aparc.$hemi.$surf.$meas
            set cmd = (mri_glmfit --table $outfile --fsgd $fsgd --glmdir $glmdir)
            echo "#C# $cmd"  | tee -a $LF
            $cmd | tee -a $LF
            if($status) exit 1;
          endif
        else
          set cmd = (asegstats2table --subjects $slist -m mean -t $outfile \
            --no-segno 0 1000 2000 --statsfile=$hemi.w-g.pct.stats $Replace53)
          echo "#C# $cmd"  | tee -a $LF
          $cmd | tee -a $LF
          if($status) exit 1;
          set glmdir = $outdir/rois/glm.$aparc.$hemi.$surf.$meas
          set cmd = (mri_glmfit --table $outfile --fsgd $fsgd --glmdir $glmdir)
          echo "#C# $cmd"  | tee -a $LF
          $cmd | tee -a $LF
          if($status) exit 1;
        endif
      end # meas
    end # hemi
  end # surf
end # aparc

endif #DoAparcStats

# ASEG Stats ----------------------------------------------------
if($DoAsegStats) then
foreach hemi ($hemilist other)
  set segidlist = `fs_temp_file --suffix .dat`
  if($hemi == lh) then
    grep Left  $ASegLUT | grep -v ypo | grep -v \# | \
    grep -v Cerebral-White-Matter | grep -v Cerebral-Cortex |\
      awk '{if($1 != 0 && length($1) !=0) print $1}' > $segidlist
  endif
  if($hemi == rh) then
    grep Right $ASegLUT | grep -v ypo | grep -v \# | \
      grep -v Cerebral-White-Matter | grep -v Cerebral-Cortex |\
      awk '{if($1 != 0 && length($1) !=0) print $1}' > $segidlist
  endif
  if($hemi == other) then 
    grep -v Left $ASegLUT | grep -v Right | \
      grep -v non | grep -v \# | \
      awk '{if($1 != 0 && length($1) !=0) print $1}' > $segidlist
  endif
  foreach meas (volume mean)
    set outfile =  $outdir/rois/aseg.$hemi.$meas.dat
    if($meas == mean) set outfile = $outdir/rois/aseg.$hemi.intensity.dat
    if(! -e $outfile || $OverWrite) then
      set cmd = (asegstats2table --subjects $slist -m $meas -t $outfile \
        --segno `cat $segidlist` $Replace53)
      echo "#C# $cmd" | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1
      set glmdir = $outdir/rois/glm.aseg.$hemi.$meas
      set cmd = (mri_glmfit --table $outfile --fsgd $fsgd --glmdir $glmdir)
      echo "#C# $cmd"  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
    endif
  end # meas
  rm -f $segidlist
end # hemi
endif 

# WMParc Stats --------------------------------------------
if($DoWMParcStats) then
set segids = (`grep -v \# $WMParcLUT | awk '{if($1 >= 3001 && $1 <= 3035 && $1 != 3004) print $1; if($1 >= 4001 && $1 <= 4035 && $1 != 4004) print $1;}'`)
set outfile =  $outdir/rois/wmparc.vol.dat
if(! -e $outfile || $OverWrite) then
  set cmd = (asegstats2table --subjects $slist -m volume -t $outfile \
    --segno $segids --stats=wmparc.stats $Replace53)
  echo "#C# $cmd" | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1
  set glmdir = $outdir/rois/glm.wmparc.vol
  set cmd = (mri_glmfit --table $outfile --fsgd $fsgd --glmdir $glmdir)
  echo "#C# $cmd"  | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif
endif

# Maps and GLM ---------------------------------------
if($DoMaps) then
  foreach meas ($mapmeaslist)
    set fmt = ()
    if($meas == w-g.pct) then
      set fmt = .mgh
      if(-e $SUBJECTS_DIR/$slist[1]/surf/$hemilist[1].w-g.pct.mgz)    set fmt = .mgz
      if(-e $SUBJECTS_DIR/$slist[1]/surf/$hemilist[1].w-g.pct.nii.gz) set fmt = .nii.gz
    endif
    foreach hemi ($hemilist)
      # concat data ---------------------
      set outfile =  $outdir/maps/$meas.$hemi.sm00.mgh
      if(! -e $outfile || $OverWrite) then
        set cmd = (mris_preproc --meas $meas$fmt --hemi $hemi --o $outfile \
           --f $outdir/subjectslist.txt --srcsurfreg $SrcSurfReg)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
      endif
      foreach fwhm ($fwhmlist)
        # smooth it -----------------
        set fwhmstr = `printf %02d $fwhm`
        set smfile =  $outdir/maps/$meas.$hemi.sm$fwhmstr.mgh
        if(! -e $smfile || $OverWrite) then
          set cmd = (mris_fwhm --no-detrend --fwhm $fwhm --i $outfile --o $smfile \
             --cortex --s fsaverage --hemi $hemi --smooth-only --prune)
          echo "#C# $cmd"  | tee -a $LF
          $cmd | tee -a $LF
          if($status) exit 1;
        endif
        # glm ------------------------
        set glmdir = $outdir/maps/glm.$meas.$hemi.sm$fwhmstr
        set cmd = (mri_glmfit --surface fsaverage $hemi --y $smfile --fsgd $fsgd --glmdir $glmdir --eres-save)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
	if($fwhm != 0) then
          set cmd = (mri_glmfit-sim --glmdir $glmdir --cwp .05 --cache 1.3 abs --2spaces)
          echo "#C# $cmd"  | tee -a $LF
          $cmd | tee -a $LF
          #if($status) exit 1;
        endif
      end # fwhm
    end # hemi
  end #meas
endif

#========================================================
# Done

echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "groupstats-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "groupstats Done" |& tee -a $LF

exit 0

###############################################

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

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

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

    case "--sd":
      if($#argv < 1) goto arg1err;
      setenv SUBJECTS_DIR $argv[1]; shift;
      breaksw

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

    case "--f":
      if($#argv < 1) goto arg1err;
      set slist = `cat $argv[1]`; shift;
      breaksw

    case "--fsgd":
      if($#argv < 1) goto arg1err;
      set fsgd = $argv[1]; shift;
      if(! -e $fsgd) then
        echo "ERROR: cannot find $fsgd"
        exit 1;
      endif
      set n = `cat $fsgd | awk '{print $1}' | grep -i contrast | wc -l`
      if($n == 0) then
        echo "ERROR: $fsgd does not have a contrast line"
        exit 1;
      endif
      breaksw

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

    case "--no-maps":
      set DoMaps = 0;
      set fwhmlist = 0;
      breaksw

    case "--lh":
      set hemilist = (lh)
      breaksw

    case "--rh":
      set hemilist = (rh)
      breaksw

    case "--new":
      set mapmeaslist = (thickness.new.mris_make_surfaces area.new.mris_make_surfaces volume)
      breaksw

    case "--base":
      # base does not have w-g.pct
      set mapmeaslist = (thickness area volume curv sulc); #  exclude w-g.pct
      breaksw

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

    case "--no-aparcstats":
      set DoAparcStats = 0
      breaksw

    case "--no-asegstats":
      set DoAsegStats = 0
      breaksw

    case "--no-wmparcstats":
      set DoWMParcStats = 0
      breaksw

    case "--no-stats":
      set DoAparcStats = 0
      set DoAsegStats = 0
      set DoWMParcStats = 0
      breaksw

    case "--keep53":
    case "--no-replace53":
      set Replace53 = ""
      breaksw

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

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

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      breaksw

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

    case "--nocleanup":
      set cleanup = 0;
      breaksw

    case "--cleanup":
      set cleanup = 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($#fsgd && $#slist) then
#  echo "ERROR: cannot spec both --fsgd and --f"
#  exit 1;
# endif
if($#fsgd == 0 && $#slist == 0) then
  echo "ERROR: must spec --fsgd or --f"
  exit 1;
endif
if($#slist == 0) set slist = `grep Input $fsgd | awk '{print $2}'`
if($#hemilist == 0) set hemilist = (lh rh)
if($#outdir == 0) then
  echo "ERROR: Must specify output folder"
  exit 1;
endif

if($#fwhmlist == 0) then
  echo "ERROR: you must specify a fwhm"
  echo "If you do not want to smooth, use --fwhm 0"
  exit 1;
endif

foreach s ($slist)
  set subjdir = $SUBJECTS_DIR/$s
  if(! -e $subjdir) then
    if($Wait) continue
    echo "ERROR: cannot find $subjdir"
    exit 1;
  endif
  if(-e $subjdir/scripts/recon-all.error) then
    echo "ERROR: $s has an error file"
    exit 1;
  endif
  set done = $subjdir/scripts/recon-all.done
  set IsRunningLH   = $subjdir/scripts/IsRunning.lh
  set IsRunningRH   = $subjdir/scripts/IsRunning.rh
  set IsRunningLHRH = $subjdir/scripts/IsRunning.lh+rh
  if(-e $done && (-e $IsRunningLH || -e $IsRunningRH || -e $IsRunningLHRH)) then
    echo "ERROR: $s has both done and IsRunning files"
    exit 1;
  endif
  if(! $Wait) then
    if(! -e $done && !(-e $IsRunningLH || -e $IsRunningRH || -e $IsRunningLHRH) ) then
      echo "ERROR: $s has neither done nor IsRunning files"
      exit 1;
    endif
  endif
end

if($Wait) then
  while (1)
    @ nwait = 0
    foreach s ($slist)
      set subjdir = $SUBJECTS_DIR/$s
      set done = $subjdir/scripts/recon-all.done
      if(! -e $done) @ nwait = $nwait + 1
    end
    if($nwait == 0) break;
    echo "Waiting for $nwait cases to finish `date`"
    sleep $tpollsec
  end
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
############--------------##################

############--------------##################
usage_exit:
  echo ""
  echo "groupstats --help"
  echo "  --o outdir : output folder"
  echo "  --fsgd group.fsgd   OR   --f subjectfile"
  echo "     if both are specified, gets subjects from subjectsfile"
  echo "  --fwhm fwhm <--fwhm fwhm2> : specify smoothing level(s)"
  echo "  --no-maps : only analyze ROI data"
  echo "  --lh : only analyze lh"
  echo "  --rh : only analyze rh"
  echo "  --sd SUBJECTS_DIR"
  echo "  --m map : use given mapname "
  echo "  --srcsurfreg SrcSurfReg (default is $SrcSurfReg)"
  echo "  --no-aparcstats : do not do aparcstats"
  echo "  --no-asegstats : do not do asegstats"
  echo "  --no-wparcstats : do not do wmparcstats"
  echo "  --no-stats : do not do any ROI stats"
  echo "  --new : append .new.mris_make_surfaces to map names"
  echo "  --base : sets meas thickness area volume curv sulc (excludes w-g.pct)" 
  echo "  --keep53 : keep 5.3 aseg names (eg, Thalamus-Proper)"
  echo "  --wait tpollsec: wait until all subjects have a .done file (poll every tpoll sec)"
  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 is a script that runs a comprehensive group analysis on both maps
and ROI results.  This script was originally designed to assist in the
testing of FreeSurfer output (in conjunction with groupstatsdiff), so
there are a lot of tests performed, more than one would normally do
for a typical study. However, it is possible to use this to do almost
all of your group analysis as it runs mris_preproc for the maps and
asegstats2table/aparcstats2table for the ROIs and then mri_glmfit for
both.

Specify the subjects through an FSGD file (--fsgd) or a subject list
file (--f) but not both. If a subject list file is specified, a
one-sample-group-mean (OSGM) test will be performed. If an FSGD file
is passed, it must have a "Contrast" line. Eg, if there are two
groups, AD and HC, then add "Contrast ad-hc 1 -1" will compute the
contrast between them. This is an alternative to using a contrast
matrix in mri_glmfit.

You can choose multiple smoothing levels with multiple --fwhm
arguments.

If the intent is to eventually use groupstatsdiff to test for a
difference in FS versions (or some other difference), then make sure
that the subject names are the same across versions.