1#!/bin/tcsh 2 3 4# - basis functions should be consistent across classes 5# i.e. should not mix GAM/BLOCK 6# - SPMG1 is okay, but no multi-component functions (AM/DM/IM are okay) 7 8# generate PPI regressors: no 3dD commands 9 10# inputs 11# - stim timing files/labels/durations 12# - PPI label (e.g. Lamy.test1) 13# - seed time series (afni_proc.py can generate example) 14# - NT per run, TR, TRnup 15 16set stim_files = ( stimuli/AV*.txt ) 17set stim_labs = ( vis aud ) 18set stim_dur = ( 20 20 ) 19 20set seed = ppi.seed.1D 21set plabel = Laud 22 23set basis = BLOCK 24 25set NT = ( 150 150 150 ) 26set TR = 2.0 27set TRnup = 20 28 29# compute some items 30# rcr - validate TRup (TR must be an integral multiple of TRup) 31set TRup = 0.1 32set demean_psych = 0 # usually 0 (for comparison, should not matter) 33 34set nstim = $#stim_files 35set run_lens = ( 300 300 300 ) 36 37set workdir = work.$plabel 38set timingdir = timing.files 39 40# ================================================================= 41# create work directory, copy inputs, enter 42 43if ( -d $workdir ) then 44 45echo "** removing old work dir $workdir..." 46rm -fr $workdir 47 # rcr - put these back 48 # echo "** will not overwrite PPI work directory $workdir, failing..." 49 # exit 1 50endif 51 52# 53mkdir $workdir 54mkdir $workdir/$timingdir 55 56cp -pv $stim_files $workdir/$timingdir 57cp -pv $seed $workdir 58set seed = $seed:t 59 60set bind = 0 61cd $workdir 62 63 64# ================================================================= 65# generate HRF 66 67# rcr - compute NT based on TRup 68if ( $basis == GAM ) then 69 70 set hrf_file = x.GAM.1D 71 3dDeconvolve -nodata 120 0.1 -polort -1 \ 72 -num_stimts 1 \ 73 -stim_times 1 1D:0 GAM \ 74 -x1D $hrf_file -x1D_stop 75 76else if ( $basis == BLOCK ) then 77 78 set hrf_file = x.BLOCK.1D 79 3dDeconvolve -nodata 150 0.1 -polort -1 \ 80 -num_stimts 1 \ 81 -stim_times 1 1D:0 "BLOCK(0.1,1)" \ 82 -x1D $hrf_file -x1D_stop 83 84else 85 echo "** invalid basis $basis, should be BLOCK or GAM (or SPMG1)" 86 exit 1 87endif 88 89 90# ================================================================= 91# create timing partition files 92 93@ bind ++ 94set prefix = p$bind.$plabel 95set timing_prefix = $prefix 96 97foreach sind ( `count -digits 1 1 $nstim` ) 98 set sind2 = `ccalc -form '%02d' $sind` 99 set tfile = $timingdir/$stim_files[$sind]:t 100 set label = $stim_labs[$sind] 101 102 if ( ! -f $tfile ) then 103 echo "** missing timing file $tfile" 104 exit 1 105 endif 106 107 timing_tool.py -timing $tfile \ 108 -tr $TRup -stim_dur $stim_dur[$sind] \ 109 -run_len $run_lens \ 110 -min_frac 0.3 \ 111 -timing_to_1D $timing_prefix.$sind2.$label \ 112 -per_run_file -show_timing 113 114 # optionally replace psychological variables with de-meaned versions 115 if ( $demean_psych ) then 116 set mean = `cat $timing_prefix.$sind2.* | 3dTstat -prefix - 1D:stdin\'` 117 echo "-- mean of psych '$label' files = $mean" 118 foreach file ( $timing_prefix.$sind2.$label* ) 119 1deval -a $file -expr "a-$mean" > rm.1D 120 mv rm.1D $file 121 end 122 endif 123end 124 125 126# ================================================================= 127# upsample seed 128 129@ bind ++ 130set prefix = p$bind.$plabel 131 132# break into n runs 133 134@ rstart = -$NT[1] 135@ rend = - 1 136foreach rind ( `count -digits 1 1 $#NT` ) 137 @ rstart += $NT[$rind] 138 @ rend += $NT[$rind] 139 1dcat $seed"{$rstart..$rend}" | 1dUpsample $TRnup stdin: \ 140 > $prefix.seed.$TRnup.r$rind.1D 141end 142 143set seed_up = $prefix.seed.$TRnup.rall.1D 144cat $prefix.seed.$TRnup.r[0-9]*.1D > $seed_up 145 146# ================================================================= 147# deconvolve 148 149set pprev = $prefix 150@ bind ++ 151set prefix = p$bind.$plabel 152set neuro_prefix = $prefix 153 154foreach rind ( `count -digits 1 1 $#NT` ) 155 3dTfitter -RHS $pprev.seed.$TRnup.r$rind.1D \ 156 -FALTUNG $hrf_file temp.1D 012 -2 \ 157 -l2lasso -6 158 1dtranspose temp.1D > $prefix.neuro.r$rind.1D 159end 160 161 162# =========================================================================== 163# partition neuro seeds 164 165set pprev = $prefix 166@ bind ++ 167set prefix = p$bind.$plabel 168 169foreach sind ( `count -digits 1 1 $nstim` ) 170 set sind2 = `ccalc -form '%02d' $sind` 171 set slab = $sind2.$stim_labs[$sind] 172 173 foreach rind ( `count -digits 1 1 $#NT` ) 174 set neuro_seed = $neuro_prefix.neuro.r$rind.1D 175 set rind2 = `ccalc -form '%02d' $rind` 176 @ nt = $NT[$rind] * $TRnup 177 178 # note partition files: 1 input, 2 outputs 179 set stim_part = $timing_prefix.${slab}_r$rind2.1D 180 set neuro_part = $prefix.a.$slab.r$rind.neuro_part.1D 181 set recon_part = $prefix.b.$slab.r$rind.reBOLD.1D 182 183 1deval -a $neuro_seed -b $stim_part -expr 'a*b' > $neuro_part 184 185 waver -FILE $TRup $hrf_file -input $neuro_part -numout $nt > $recon_part 186 end 187 188 # and generate upsampled seeds that span runs 189 cat $prefix.b.$slab.r*.reBOLD.1D > $prefix.d.$slab.rall.reBOLD.1D 190end 191 192# and generate corresponding (reBOLD) seed time series 193foreach rind ( `count -digits 1 1 $#NT` ) 194 set neuro_seed = $neuro_prefix.neuro.r$rind.1D 195 waver -FILE $TRup $hrf_file -input $neuro_seed -numout $nt \ 196 > $prefix.c.seed.r$rind.reBOLD.1D 197end 198 199# to compare with $seed_up 2003dMean -sum -prefix - $prefix.d.[0-9]*.1D > $prefix.d.task.rall.reBOLD.1D 201cat $prefix.c.seed.r*.reBOLD.1D > $prefix.d.seed.rall.reBOLD.1D 202echo == can compare upsampled seeds: \ 203 $seed_up $prefix.d.{seed,task}.rall.reBOLD.1D 204set seed_rebold_up = $prefix.d.seed.rall.reBOLD.1D 205 206 207# =========================================================================== 208# downsample 209 210set pprev = $prefix 211@ bind ++ 212set prefix = p$bind.$plabel 213 214foreach rind ( `count -digits 1 1 $#NT` ) 215 set neuro_seed = $neuro_prefix.neuro.r$rind.1D 216 @ nt = $NT[$rind] * $TRnup 217 218 foreach sind ( `count -digits 1 1 $nstim` ) 219 set sind2 = `ccalc -form '%02d' $sind` 220 set recon_part = $pprev.b.$sind2.$stim_labs[$sind].r$rind.reBOLD.1D 221 set recon_down = $prefix.$sind2.$stim_labs[$sind].r$rind.PPIdown.1D 222 223 1dcat $recon_part'{0..$('$TRnup')}' > $recon_down 224 end 225 226 # and downsample filtered seed time series 227 1dcat $seed_rebold_up'{0..$('$TRnup')}' > $seed:r.reBOLD.1D 228end 229 230 231# =========================================================================== 232# catentate across runs: final PPI regressors 233 234set pprev = $prefix 235@ bind ++ 236set prefix = p$bind.$plabel 237 238foreach sind ( `count -digits 1 1 $nstim` ) 239 set sind2 = `ccalc -form '%02d' $sind` 240 set slab = $sind2.$stim_labs[$sind] 241 242 cat $pprev.$slab.r*.PPIdown.1D > $prefix.$slab.rall.PPI.1D 243end 244 245# ================================================================= 246# make a final comparison time series 247 248set pprev = $prefix 249@ bind ++ 250set prefix = p$bind.$plabel 251 2523dMean -sum -prefix - $pprev.* > $prefix.sum.PPI.1D 253 254echo "== can compare original seed to sum of PPI regressors:" 255echo " 1dplot -one $seed $prefix.sum.PPI.1D" 256 257echo "" 258echo "== final PPI regressors: " $seed $pprev.* 259echo " (copy to stimuli dir)" 260echo "" 261