#!/bin/tcsh

# ----------------------------------------------------------------------
# compute GCOR (average global correlation)
#
# GCOR = average correlation of every voxel with every other voxel
#        (i.e. average of the correlation matrix)
#
#    inputs: an EPI dataset and a mask (on same grid, of course)
#    output: a single GCOR value
# ----------------------------------------------------------------------


# ----------------------------------------------------------------------
# main variables

# dataset inputs
set inset       = ''            # input dataset
set mset        = ''            # mask dataset
set corrset     = ''            # correlation dataset

set demean      = 1             # whether to demean data (faster if not)
set nfirst      = 0             # number of first TRs to remove
set savetmp     = 0             # save temporary files?
set tmpunit     = tmp.unit      # 3dTnorm output
set tmpave      = tmp.unit.gmean.1D # 3dmaskave output
set tmpgc       = tmp.gcor.1D   # GCOR result (stupid stderr)
set tmpprod     = tmp.uprod     # unit time series times gmean time series
set verb        = 1             # verbose level (0 = quiet)

set prog = @compute_gcor

set version = "0.0, 17 Jan, 2013"

if ( $#argv < 1 ) goto SHOW_HELP

# ----------------------------------------------------------------------
# check for options
set ac = 1
while ( $ac <= $#argv )
   if      ( "$argv[$ac]" == "-help" ) then
      goto SHOW_HELP
   else if ( "$argv[$ac]" == "-hist" ) then
      goto SHOW_HIST
   else if ( "$argv[$ac]" == "-ver" ) then
      echo "version $version"
      exit

   else if ( "$argv[$ac]" == "-corr_vol" ) then
      @ ac += 1
      if ( $#argv < $ac ) then
         echo "** missing argument for -corr_vol"
         exit 1
      endif
      set corrset = $argv[$ac]
   else if ( "$argv[$ac]" == "-input" ) then
      @ ac += 1
      if ( $#argv < $ac ) then
         echo "** missing argument for -input"
         exit 1
      endif
      set inset = $argv[$ac]
   else if ( "$argv[$ac]" == "-mask" ) then
      @ ac += 1
      if ( $#argv < $ac ) then
         echo "** missing argument for -mask"
         exit 1
      endif
      set mset = $argv[$ac]
   else if ( "$argv[$ac]" == "-nfirst" ) then
      @ ac += 1
      if ( $#argv < $ac ) then
         echo "** missing argument for -nfirst"
         exit 1
      endif
      set nfirst = $argv[$ac]
   else if ( "$argv[$ac]" == "-no_demean" ) then
      set demean = 0
   else if ( "$argv[$ac]" == "-savetmp" ) then
      set savetmp = 1
   else if ( "$argv[$ac]" == "-verb" ) then
      @ ac += 1
      if ( $#argv < $ac ) then
         echo "** missing argument for -verb"
         exit 1
      endif
      set verb = $argv[$ac]
      if ( $verb > 2 ) then
         echo "++ turning echo on..."
         set echo
      endif
   else
      echo "** unknown option: $argv[$ac]"
      exit 1
   endif
   @ ac += 1
end

# ----------------------------------------------------------------------
# check for inputs
if ( $inset == '' ) then
   echo "** missing input dataset, see -help output for details"
   exit 1
endif

# quick test on input dataset (and possibly mask), while setting view
set view = `3dinfo -av_space $inset`
if ( $status || $view == NO-DSET ) then
   echo "** failed to get view of -input $inset, check command"
   exit 1
endif

if ( $mset != '' ) then
   set vv = `3dinfo -av_space $inset`
   if ( $status || $vv == NO-DSET ) then
      echo "** failed to get view of -mask $mset, check command"
      exit 1
   endif
endif

# ----------------------------------------------------------------------
# get to work
# ----------------------------------------------------------------------

# demean?
set polopt = ''
if ( $demean ) set polopt = '-polort 0'

# ----------------------------------------------------------------------
# if 1D input, just run 1d_tool.py
if ( $inset:e == 1D || $inset:e == txt ) then
   if ( $mset != '' ) then
      echo "** -mask is not currently allowed with 1D input (pester Rick)"
      exit 1
   endif
   if ( $nfirst > 0 ) then
      echo "** -nfirst is not currently allowed with 1D input (pester Rick)"
      exit 1
   endif
   if ( $verb > 1 ) echo 1d_tool.py -verb $verb -infile $inset -show_gcor
   1d_tool.py -verb $verb -infile $inset -show_gcor
   exit $status
endif

# ----------------------------------------------------------------------
# create unit dataset

# can separate command, but nfirst might need care...
set cmd = ( 3dTnorm -overwrite $polopt -prefix $tmpunit )
if ( $verb > 0 ) then
   echo $cmd $inset\[$nfirst..\$]
   $cmd $inset\[$nfirst..\$]
else
   ( $cmd $inset\[$nfirst..\$] ) >& /dev/null
endif

if ( ! -f $tmpunit$view.HEAD ) then
   echo "** command failure:" $cmd $inset\[$nfirst..\$]
   exit
endif

# ----------------------------------------------------------------------
# get global mean of unit time series
set mopt = ''
if ( $mset != "" ) set mopt = ( -mask $mset )

# this time, take care with redirection
set cmd = ( 3dmaskave -quiet $mopt $tmpunit$view )
if ( $verb > 0 ) then
   echo $cmd \> $tmpave
   $cmd > $tmpave
else
   ( $cmd > $tmpave ) >& /dev/null
endif

# ----------------------------------------------------------------------
# if requested, compute the correlation volume
if ( $corrset != "" ) then
   # messy to include quotes or * this way...
   set c1 = ( 3dcalc -a $tmpunit$view -b $tmpave )
   set c2 = ( 3dTstat -sum -prefix $corrset $tmpprod$view )
   if ( $verb > 0 ) then
      echo $c1 -expr a\*b -prefix $tmpprod
           $c1 -expr a\*b -prefix $tmpprod
      echo $c2
           $c2
   else
      $c1 -expr a\*b -prefix $tmpprod >& /dev/null
      $c2 >& /dev/null
   endif
endif
   
# ----------------------------------------------------------------------
# and finally GCOR, the squared length of that average unit time series
# (which is not itself, unit)

set cmd = ( 3dTstat -sos -prefix - )
if ( $verb > 0 ) then
   echo $cmd $tmpave\'
   $cmd $tmpave\' > $tmpgc
   echo GCOR = `cat $tmpgc`
else
   ( $cmd $tmpave\' > $tmpgc ) >& /dev/null
   cat $tmpgc
endif
   
   
# ----------------------------------------------------------------------
# cleanup
if ( ! $savetmp ) then
   if ( $verb > 1 ) echo rm -f $tmpunit$view\* $tmpave $tmpgc
   \rm -f $tmpunit$view* $tmpave $tmpgc 
   if ( -f $tmpprod$view.HEAD ) then
      \rm -f $tmpprod$view.*
   endif
else if ( $verb > 1 ) then
   echo "keeping temp files: $tmpunit$view* $tmpave $tmpgc"
   if ( -f $tmpprod$view.HEAD ) then
      echo "keeping temp files: $tmpprod$view.*"
   endif
endif

# ======================================================================
# DONE
# ======================================================================

exit


SHOW_HELP:
# ----------------------------------------------------------------------
cat << EOF

-----------------------------------------------------------------
$prog       - compute GCOR, the global correlation

    usage : $prog [options] -input dataset

This program computes the average correlation between every voxel
and ever other voxel, over any given mask.  This output GCOR value
is a single number.

-----------------------------------------------------------------

Common examples:

  0. This program can be used for 1D files:

        $prog -input data.1D

     HOWEVER, if column selection is desired, please use 1d_tool.py, directly.

        1d_tool.py -infile data.1D'[2..17]' -show_gcor

  1. Simple usage, akin to the afni_proc.py processing script.

        $prog -input errts.FT+orig -mask full_mask.FT+orig

     OR, for +tlrc:

        $prog -input errts.FT+tlrc -mask full_mask.FT+tlrc

  2. Speed things up slightly, an errts dataset does not need to be demeaned.

        $prog -no_demean -input errts.FT+tlrc -mask full_mask.FT+tlrc

  3. Be vewy, veeewy, qwiet...

        $prog -verb 0 -input errts.FT+tlrc -mask full_mask.FT+tlrc

     OR, save the result:

        set gcor = \`$prog -verb 0 -input errts.FT+tlrc -mask full_mask.FT+tlrc\`

  4. Output correlation volume: for each voxel, the average correlation
     with all voxels in mask.

     Specify correlation volume prefix, FT_corr.

       $prog -input errts.FT+tlrc -mask full_mask.FT+tlrc -corr_vol FT_corr

-----------------------------------------------------------------

Overview of processing steps: 

0.  If the input is a 1D file, cheat and run "1d_tool.py -show_gcor", instead.

otherwise...

1.  Scale the input to a unit time series, so that each voxel voxel has a
    length of 1.

2.  Compute the average of these unit time series.

3a. If requested, compute the correlation volume, the dot product of the
    unit and average time series.

3b.  Return GCOR = the length of the resulting average, squared.

---------------------------------------------

terminal options:

   -help             : show this help
   -hist             : show modification history
   -ver              : show version number

important processing options:

   -input DSET       : specify input dataset to compute the GCOR over
   -mask DSET        : specify mask dataset, for restricting the computation

other processing options:

   -corr_vol PREFIX  : specify input dataset to compute the GCOR over
   -nfirst NFIRST    : specify number of initial TRs to ignore
   -no_demean        : do not (need to) demean as first step
   -savetmp          : save temporary files (do not remove at end)
   -verb VERB        : set verbose level (0=quiet, 3=max)

---------------------------------------------

R Reynolds, Jan, 2013
------------------------------------------------------------
EOF

exit

SHOW_HIST:

cat << EOF
-----------------------------------------------------------------
$prog modification history:

   0.0  : Jan 17, 2013: initial version
   0.1  : Feb 27, 2015: added -corr_vol

EOF

exit