#!/bin/tcsh

# alphasurface.tcsh

# This script runs a simulation to find the cluster correction size in order to correct for multiple comparisons for surface data. It populates the input surface with random data, applies spatial smoothing to those data to a target FWHM, clusters the data, and computes the maximum in area or number of surface vertices (change the variable as desired). The simulation is designed to correct for the specified Family Wise Error rate.
  
# By Matthew Schiel, Michael Andric, and Anthony Dick; 4/23/10.

# Node = Vertex = Surface vertex

# This script assumes you are running from the SUMA directory with the spec file and surfaces.

set iter = "1000" # number of iterations
set fwhm = "2" # smooth to the target FWHM; find target FWHM from actual data using SurfFWHM 
set pthr = ".001" # set per vertex threshold
set rmm  = "-2" # see help file for SurfClust
set spec = "010_lh.spec" # the spec file
set surf = "pial" # the surface file
set Nnode = "A" # node automatic; this will work unless the SUMA file formats change. Double check the accuracy in the screen output.
set sortby = "area" #choose "n_nodes" to simulate maximum cluster size by number of nodes; choose "area" to simulate cluster size by area 
set FWEr = ".05" # Family-Wise Error rate correction alpha 
set prefix = "testingsimulation" # the output filename

### Don't change anything below this line unless you know what you are doing

if ( "$Nnode" == "A") then
    set surfpath = `grep -B 2 "SurfaceState = $surf" "$spec" | head -n 1 | awk '{print $3}'` # Find the SurfaceState text in the spec file that points to the surface listed in $surf, and print the filename two lines above that
    set Nnode = `grep -v -m 1 "^#" "$surfpath" | awk '{print $1}'` #print the number of surface vertices in that filename, skipping the commented text
    echo "Nnode is now automatically set to $Nnode"
endif

set t = 0
while ( $t < $iter )
	echo "Simulation number is `expr $t + 1` out of $iter"
    1deval -num $Nnode -del 1 -expr 'gran(0,1)' > __temp_rand.1D.dset # generate random data from a normal distribution
    SurfSmooth -spec "$spec" -surf_A "$surf" -met HEAT_07 -target_fwhm "$fwhm" -input __temp_rand.1D.dset -output __temp_rand.smooth.1D.dset
    sort -g __temp_rand.smooth.1D.dset > __temp_rand.smooth.sort.1D.dset # Smooth those data with given input parameters
    set threshindex = `ccalc -int -expr "(1-$pthr)*$Nnode - 1"` # take the percentile of the number of vertices given by 1-pthr (e.g., 1 -.01 = 99%) and adjust row numbering to meet AFNI specifications
    set threshvalue = `1dcat -a "__temp_rand.smooth.sort.1D.dset{$threshindex}"` #give the value at that percentile from the distribution of smoothed, random data
    echo "After smoothing the z-distributed simulated data, your threshold on the new distribution is $threshvalue" 
    SurfClust -spec "$spec" -surf_A "$surf" -input __temp_rand.smooth.1D.dset 0 -rmm "$rmm" -thresh_col 0 -thresh "$threshvalue" -n 1 -sort_$sortby -no_cent -prefix __temp_rand #Cluster the smoothed random data
    if ( "$t" == "0" ) then
	grep -B 1 -v -m 1 "#" __temp_rand_ClstTable* > "$prefix".1D # grab column headers and first cluster line
    else
	grep -v -m 1 "#" __temp_rand_ClstTable* >> "$prefix".1D # ignore column headers after first iteration
    endif
    rm -f __temp_rand*
    @ t = $t + 1
end
if ( "$sortby" == "area" ) then
	set sortcol = 3
else
	set sortcol = 2
endif	
grep -v "#" "$prefix".1D | sort -g -k $sortcol > "$prefix".sorted.1D # Sort the output by the sortby parameter (area or nodes), and output file without comments
set flength = `grep '.' "$prefix".sorted.1D | wc -l  | awk '{print $1}'` # find the filelength not counting blank lines
set threshindex = `ccalc -int -expr "(1-$FWEr)*$flength - 1"` # find the index of the value at FWE correction percentile, adjust row numbering to meet AFNI specifications   
@ sortcol = $sortcol - 1 # adjust to meet AFNI specifications
set threshvalue = `1dcat -a "$prefix.sorted.1D{$threshindex}[$sortcol]"` # find the value at FWE correction percentile 
echo "To control for the Family Wise Error rate of $FWEr, the minimum cluster size in $sortby is $threshvalue" # give result