#!/bin/tcsh
# SEMscript.csh
# 
# Compute structural equation model for subject, ROI data

# need minimum of at least 2 ROIs, 2 Subjects for this script

# Each subject has a single input dataset composed of the columns
#  corresponding to each ROI
# for filenames without spaces

# call this script with 
#
# tcsh SEMscript.csh subj1_allrois.1D subj2_allrois.1D subj3_allrois3.1D ... \
#   MinModel.1D
#
# Output includes Structural Equation Model coefficients for fixed paths (model
#   validation), variable number of paths starting with a minimum model
#   and grown two different ways to a maximum number of coefficients (tree and
#   forest growth), correlation, covariance matrices and corresponding 
#   degrees of freedom and the residual variance
#

# set overall parameters here
@ last = $#argv - 1

set minModel = $argv[$#argv]
set subj_dsets = ( $argv[1-$last] )

# set subj_dsets = `glob $1`
# set minModel = $2
set NSUBJS = $#subj_dsets
glob $subj_dsets ; echo
# number of time points is from number of elements (lines) in each column
set NVALS = `1dcat ${subj_dsets[1]}'[0]' | wc -l`
# number of ROIs is number of elements (words) in the first (and following) row
set NROIS = `1dcat ${subj_dsets[1]}'{0}' | wc -w`

@ NP = $NVALS + 2
@ NVM1 = $NVALS - 2

# remove previous files
rm -f SEM_roi*_extract.1D
rm -f SEM_roi*_allsubjs.1D
rm -f SEM_svd_roi*.txt
rm -f SEM_eigents_roi*.1D
rm -f SEM_sv_roi*.1D
rm -f SEM_sv2_roi*.1D
rm -f SEM_resvar.1D
rm -f SEM_resvarcorr.1D
rm -f SEM_resvar_roi*.1D
rm -f SEM_ts_mean_roi*.1D
rm -f SEM_eigentsc_roi*.1D
rm -f SEM_autocorr.1D
rm -f SEM_cov.1D
rm -f SEM_corr.1D
rm -f SEM_temp.1D
rm -f SEM_DF.1D

foreach roi (`count 1 $NROIS`)
   @ roi_index = $roi - 1
   foreach subj (`count 1 $NSUBJS`)
      1dcat ${subj_dsets[${subj}]}"[${roi_index}]" >  \
        SEM_roi${roi}_subj${subj}_extract.1D
   end
   1dcat SEM_roi${roi}_subj*_extract.1D > SEM_roi${roi}_allsubjs.1D

   # singular value decomposition for eigentimeseries and singular values
   
   1dsvd SEM_roi${roi}_allsubjs.1D > SEM_svd_roi${roi}.txt
   grep -A $NP -m 1 "++ Left Vectors" SEM_svd_roi${roi}.txt | \
     tail -$NVALS | tr -d ":" > SEM_temp.1D
   1dcat SEM_temp.1D'[1]' > SEM_eigents_roi${roi}.1D


   #  compute residual variance   
#   grep -A 1 -m 1 "++ Left Vectors" SEM_svd_roi${roi}.txt | \
#     tail -1 > SEM_temp.1D
#   1dtranspose SEM_temp.1D > SEM_sv_roi${roi}.1D
#   1deval -a SEM_sv_roi${roi}.1D -expr 'a*a' > SEM_sv2_roi${roi}.1D
   1dcat SEM_temp.1D'[1]' > SEM_sv_roi${roi}.1D
   1deval -a SEM_sv_roi${roi}.1D -expr 'a*a' > SEM_sv2_roi${roi}.1D

   3dTstat -sum -prefix SEM_resvar_roi${roi} SEM_sv2_roi${roi}.1D'{1..$}'\'
   1dcat SEM_resvar_roi${roi}.1D >> SEM_resvar.1D 
	
   #  compute residual variance corresponding to correlation matrix     
   set s1 = `1dcat SEM_sv2_roi${roi}.1D'{0}'`
   set ssum = `1dcat SEM_resvar_roi${roi}.1D`
   ccalc -expr "1- $s1/ ($ssum+$s1)" >> SEM_resvarcorr.1D  

   # flip sign of eigentimeseries as needed
   3dMean -prefix SEM_ts_mean_roi$roi SEM_roi${roi}_subj*_extract.1D
   set dotp = `3ddot -dodot SEM_eigents_roi${roi}.1D SEM_ts_mean_roi${roi}.1D`
   1deval -a SEM_eigents_roi${roi}.1D  \
    -expr "step(-$dotp)*(-a)+step($dotp)*a" > SEM_eigentsc_roi${roi}.1D

   # compute first order autocorrelation coefficient
   3ddot -demean SEM_eigentsc_roi${roi}.1D"{0..${NVM1}}"\
                 SEM_eigentsc_roi${roi}.1D'{1..$}' >> SEM_autocorr.1D
end

# compute covariance matrix
1ddot -dem -cov -terse SEM_eigentsc_roi*.1D > SEM_cov.1D

# compute correlation matrix
1ddot -dem -terse SEM_eigentsc_roi*.1D > SEM_corr.1D

# compute degrees of freedom, DF
1deval -a SEM_autocorr.1D -expr "(${NVALS} / ${NROIS}) * (1-a)/(1+a)" > SEM_temp.1D
3dTstat -prefix SEM_DF -sum SEM_temp.1D\'
set DF = `1dcat SEM_DF.1D`

#1dSEM -theta $minModel -C SEM_cov.1D -DF $DF -psi SEM_resvar.1D -limits -10 10
#1dSEM -theta $minModel -C SEM_corr.1D -DF $DF -psi SEM_resvarcorr.1D -limits -2 2

#1dSEM -theta $minModel -C SEM_cov.1D -DF $DF -psi SEM_resvar.1D -tree_growth -limits -10 10
1dSEM -theta $minModel -C SEM_corr.1D -DF $DF -psi SEM_resvarcorr.1D -tree_growth -limits -2 2
#1dSEM -theta $minModel -C SEM_cov.1D -DF $DF -psi SEM_resvar.1D -forest_growth -leafpicker -limits -10 10
#1dSEM -theta $minModel -C SEM_corr.1D -DF $DF -psi SEM_resvarcorr.1D -forest_growth -leafpicker -limits -2 2

# clean up temporary files
rm -f SEM_temp.1D
rm -f SEM_roi*_extract.1D
rm -f SEM_roi*_allsubjs.1D
