Index: ./doc/fiasco_executables.html *** /TOIL-U4/fiasco/Fiasco5.1/doc/fiasco_executables.html Wed Feb 13 19:01:06 2002 --- ./doc/fiasco_executables.html Thu Apr 4 13:46:11 2002 *************** *** 27,32 **** --- 27,38 ----
merge_detrend.csh
combines multiple Fiasco session results for the same subject by aligning and concatenating the detrended session data. +
pooled_mean.csh +
constructs the pooled mean of a group of means with associated + counts. +
pooled_stdv.csh +
constructs the pooled standard deviation of a group of stdvs with + associated counts.
testingroup.csh
constructs T and P scores for response to a two-condition test for multiple subjects, all of whom are in the same group. *************** *** 185,190 **** --- 191,198 ----
is called by scripts to provide their help information.
sgrid
spiral gridding +
smooth_missing
linearly interpolates over missing data in time series.
smregpar *************** *** 243,249 ****
Joel Welling
! Last modified: Thu Jan 24 18:07:26 EST 2002 --- 251,257 ----
Joel Welling
! Last modified: Wed Mar 20 15:45:47 EST 2002 Index: ./src/baseline/Makefile *** /TOIL-U4/fiasco/Fiasco5.1/src/baseline/Makefile Wed Feb 13 19:01:12 2002 --- ./src/baseline/Makefile Tue Mar 5 17:30:31 2002 *************** *** 11,17 **** PKG_MAKEBINS = $(CB)/baseline PKG_LIBS = -lfmri -lmri -lpar -lbio -lacct -lmisc -lcrg -lm - include ../Makefile.pkg LIBFILES= $L/libmri.a $L/libpar.a $L/libbio.a $L/libarray.a $L/libmisc.a \ $L/libacct.a --- 11,16 ---- *************** *** 21,26 **** --- 20,27 ---- CSOURCE= baseline.c DOCFILES= baseline_help.help + include ../Makefile.pkg + $O/baseline.o: baseline.c $(CC_RULE) Index: ./src/csh/cg_brain_stats.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/cg_brain_stats.csh Wed Feb 13 19:01:16 2002 --- ./src/csh/cg_brain_stats.csh Tue Apr 2 21:00:33 2002 *************** *** 32,38 **** # the 'interesting' parts for display. # echo '#'`date` $0 ! echo '#$Id: cg_brain_stats.csh,v 1.14 2002/01/10 18:34:06 bakalj Exp $' if ( ! -d $F_STAT_MPATH) mkdir $F_STAT_MPATH if ( ! -d $F_STAT_SPATH) mkdir $F_STAT_SPATH --- 32,38 ---- # the 'interesting' parts for display. # echo '#'`date` $0 ! echo '#$Id: cg_brain_stats.csh,v 1.15 2002/04/03 02:00:39 welling Exp $' if ( ! -d $F_STAT_MPATH) mkdir $F_STAT_MPATH if ( ! -d $F_STAT_SPATH) mkdir $F_STAT_SPATH *************** *** 72,82 **** # Make a temporary file with just fit data mri_copy_chunk -input $1 -o brain_stat_tmp -c images -t images -f .dat ! # Build a map of probability from log of Bayes factor mri_copy_chunk -input $1 -o brain_stat_lnb_tmp -c lnbfac -t images -f .dat ! mri_rpn_math '$1,exp,'${F_BAYES_P0}',*,dup,1.0,+,/' -o $F_STAT_MPATH/Bayes.P \ ! brain_stat_lnb_tmp ! rm brain_stat_lnb_tmp.dat brain_stat_lnb_tmp.mri # extract baselevel statistics @ base_stde_offset = $n_estimates --- 72,85 ---- # Make a temporary file with just fit data mri_copy_chunk -input $1 -o brain_stat_tmp -c images -t images -f .dat ! # Build a map of probability from log of Bayes factor. We convert the ! # input to doubles so that the final output will be in doubles. mri_copy_chunk -input $1 -o brain_stat_lnb_tmp -c lnbfac -t images -f .dat ! mri_type_convert -double brain_stat_lnb_tmp brain_stat_lnb_dbl_tmp ! mri_rpn_math '$1,exp,'${F_BAYES_P0}',*,dup,1.0,+,/' \ ! -o $F_STAT_MPATH/Bayes.P brain_stat_lnb_dbl_tmp ! mri_destroy_dataset brain_stat_lnb_tmp ! mri_destroy_dataset brain_stat_lnb_dbl_tmp # extract baselevel statistics @ base_stde_offset = $n_estimates *************** *** 160,166 **** endif # Clean up the temp file ! rm brain_stat_tmp.dat brain_stat_tmp.mri --- 163,169 ---- endif # Clean up the temp file ! mri_destroy_dataset brain_stat_tmp Index: ./src/csh/epi.defaults.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/epi.defaults.csh Wed Feb 13 19:01:16 2002 --- ./src/csh/epi.defaults.csh Thu Mar 21 19:03:22 2002 *************** *** 27,33 **** # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '#$Id: epi.defaults.csh,v 1.31 2001/12/21 01:44:57 welling Exp $' if(! -d in)mkdir in # # --- 27,33 ---- # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '#$Id: epi.defaults.csh,v 1.33 2002/03/22 00:04:21 welling Exp $' if(! -d in)mkdir in # # *************** *** 96,101 **** --- 96,108 ---- setenv F_READER_OPTS "" # # + # epi.despike.csh + # + setenv F_DESPIKE_OUTPUT despike + setenv F_DESPIKE_PARMS spike_count.par + setenv F_DESPIKE_CUTOFF 4.0 + # + # # epi.baseline.csh # # baseline header file *************** *** 138,143 **** --- 145,163 ---- # deghost filter bandwidth setenv F_DEGHOST_BAND 8.0 # + # phase_lock.csh + # + # phase_lock header file + setenv F_PHLOCK_OUTPUT phase_lock + # unsmoothed estimates + setenv F_PHLOCK_RWPARMS phase_lock_raw.par + # smoothed estimates + setenv F_PHLOCK_PARMS phase_lock.par + # smoother bandwidth + setenv F_PHLOCK_BAND 3.0 + # algorithm: "weighted" or "center" + setenv F_PHLOCK_ALG "weighted" + # # epi.meanc.csh # # meanc header file Index: ./src/csh/false_discovery.py *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/false_discovery.py Wed Feb 13 19:01:19 2002 --- ./src/csh/false_discovery.py Mon Apr 1 17:19:20 2002 *************** *** 157,169 **** continue cutoff= (scale*i)/validpix ! # print "%d: %f vs %f"%(i,cutoff,val) if cutoff >= val: rejectNullThresh= val i= i+1 ! print rejectNullThresh # Clean up removeTmpDir(tmpdir); --- 157,170 ---- continue cutoff= (scale*i)/validpix ! # print "%d: %g vs %g"%(i,cutoff,val) if cutoff >= val: rejectNullThresh= val i= i+1 ! # We use repr() to get more significant digits in output. ! print repr(rejectNullThresh); # Clean up removeTmpDir(tmpdir); Index: ./src/csh/FIASCO *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/FIASCO Fri Mar 1 19:03:48 2002 --- ./src/csh/FIASCO Thu Apr 4 14:08:46 2002 *************** *** 45,58 **** echo '#BEGINNING NEW ANALYSIS' >>& $LOGFILE setenv F_STARTDATE "`date`" echo '#'$F_STARTDATE $0 >>& $LOGFILE ! echo '$Id: FIASCO,v 1.22 2002/03/01 21:14:43 welling Exp $'\ >>& $LOGFILE echo '#'`whoami`' '`hostname` >>& $LOGFILE echo '#'$1' Processing' >>& $LOGFILE # FIASCO_VERSION number setenv FIASCO_VERSION 5.1 ! setenv FIASCO_PATCHLVL 1 echo "#FIASCO Version "$FIASCO_VERSION" patch "$FIASCO_PATCHLVL\ >>& $LOGFILE if ($1 == 'spiral') then --- 45,58 ---- echo '#BEGINNING NEW ANALYSIS' >>& $LOGFILE setenv F_STARTDATE "`date`" echo '#'$F_STARTDATE $0 >>& $LOGFILE ! echo '$Id: FIASCO,v 1.23 2002/04/04 19:08:35 welling Exp $'\ >>& $LOGFILE echo '#'`whoami`' '`hostname` >>& $LOGFILE echo '#'$1' Processing' >>& $LOGFILE # FIASCO_VERSION number setenv FIASCO_VERSION 5.1 ! setenv FIASCO_PATCHLVL 2 echo "#FIASCO Version "$FIASCO_VERSION" patch "$FIASCO_PATCHLVL\ >>& $LOGFILE if ($1 == 'spiral') then Index: ./src/csh/groupcompare.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/groupcompare.csh Wed Feb 13 19:01:17 2002 --- ./src/csh/groupcompare.csh Thu Apr 4 13:48:37 2002 *************** *** 3,9 **** # This script does between-group comparisons. # echo '#'`date` $0 ! echo '# $Id: groupcompare.csh,v 1.15 2002/01/10 18:34:06 bakalj Exp $' # Check for help command if ( $#argv >= 1 ) then --- 3,9 ---- # This script does between-group comparisons. # echo '#'`date` $0 ! echo '# $Id: groupcompare.csh,v 1.20 2002/04/04 18:48:39 welling Exp $' # Check for help command if ( $#argv >= 1 ) then *************** *** 173,183 **** mri_rpn_math -o rawfactor3 '$1,$2,*' rawfactor1 rawfactor2 ! set factorlist = ( rawfactor1 rawfactor2 rawfactor3 ) set effectfactor = 1 set groupfactor = 2 set crossfactor = 3 set which_subj = 0 foreach subj ( $group0_subjects ) if ( $subj == $lastgrp0 ) then --- 173,184 ---- mri_rpn_math -o rawfactor3 '$1,$2,*' rawfactor1 rawfactor2 ! set rawfactorlist = ( rawfactor1 rawfactor2 rawfactor3 ) set effectfactor = 1 set groupfactor = 2 set crossfactor = 3 + set subjfactorlist = ( ) set which_subj = 0 foreach subj ( $group0_subjects ) if ( $subj == $lastgrp0 ) then *************** *** 187,193 **** mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,$2,1.0,$3,/,*,-,$4,*' \ conds_x_subjects group0_mask group0_count rawfactor1 ! set factorlist = ( $factorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end foreach subj ( $group1_subjects ) --- 188,194 ---- mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,$2,1.0,$3,/,*,-,$4,*' \ conds_x_subjects group0_mask group0_count rawfactor1 ! set subjfactorlist = ( $subjfactorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end foreach subj ( $group1_subjects ) *************** *** 198,220 **** mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,$2,1.0,$3,/,*,-,$4,*' \ conds_x_subjects group1_mask group1_count rawfactor1 ! set factorlist = ( $factorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end ! @ nfactors = ${#factorlist} ! echo "nfactors is " $nfactors ! echo "unscaled factors are " $factorlist echo "The factor of interest is number " $crossfactor # Remap the factors into the order needed for GLM ! foreach fname ( $factorlist ) mri_remap -file $fname -order t -extents ${tdim} end ! # Make the factor part of the operand string for mri_glm ! set factoropts = ( ) ! foreach fname ( $factorlist ) ! set factoropts = ( $factoropts -f ${fname}:images ) end # Build the scaling matrix and the vector of mean values. There are --- 199,229 ---- mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,$2,1.0,$3,/,*,-,$4,*' \ conds_x_subjects group1_mask group1_count rawfactor1 ! set subjfactorlist = ( $subjfactorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end ! @ nrawfactors = ${#rawfactorlist} ! @ nsubjfactors = ${#subjfactorlist} ! @ nfactors = $nrawfactors + $nsubjfactors ! echo "number of fixed effect factors is " $nrawfactors ! echo "number of factors associated with individual subjects is " $nsubjfactors ! echo "total nfactors is " $nfactors ! echo "unscaled factors are " $rawfactorlist $subjfactorlist echo "The factor of interest is number " $crossfactor # Remap the factors into the order needed for GLM ! foreach fname ( $rawfactorlist $subjfactorlist ) mri_remap -file $fname -order t -extents ${tdim} end ! # Make the factor part of the operand strings for mri_glm ! set rawfactoropts = ( ) ! foreach fname ( $rawfactorlist ) ! set rawfactoropts = ( $rawfactoropts -f ${fname}:images ) ! end ! set subjfactoropts = ( ) ! foreach fname ( $subjfactorlist ) ! set subjfactoropts = ( $subjfactoropts -f ${fname}:images ) end # Build the scaling matrix and the vector of mean values. There are *************** *** 235,250 **** set stdv_fname = `map_name.py -d cond=$cond -d file=Stdv $subj` set counts_fname = `map_name.py -d cond=$cond -d file=Count $subj` set tfile1 = junk2_${subj}_${cond} ! # mri_rpn_math -o $tfile1 \ ! # '0.0,$1,dup,1.0,<=,if_keep,sqrt,1.0,$2,dup,0.0,!=,if_keep,/,0.0,$2,0.0,==,if_keep' \ ! # $counts_fname $stdv_fname mri_rpn_math -o $tfile1 '0.0,$1,dup,1.0,<=,if_keep,sqrt' \ $counts_fname $stdv_fname mri_remap -file $tfile1 -order vtxyz mri_delete_chunk -dataset $tfile1 -chunk missing set scale_args = ( $scale_args $tfile1 ) ! set tfile2 = junk3_${subj}_${cond} ! mri_copy_dataset ${mean_fname} ${tfile2} mri_remap -file $tfile2 -order vtxyz mri_delete_chunk -dataset $tfile2 -chunk missing set rawy_args = ( $rawy_args $tfile2 ) --- 244,256 ---- set stdv_fname = `map_name.py -d cond=$cond -d file=Stdv $subj` set counts_fname = `map_name.py -d cond=$cond -d file=Count $subj` set tfile1 = junk2_${subj}_${cond} ! set tfile2 = junk3_${subj}_${cond} mri_rpn_math -o $tfile1 '0.0,$1,dup,1.0,<=,if_keep,sqrt' \ $counts_fname $stdv_fname mri_remap -file $tfile1 -order vtxyz mri_delete_chunk -dataset $tfile1 -chunk missing set scale_args = ( $scale_args $tfile1 ) ! mri_copy_dataset ${mean_fname} $tfile2 mri_remap -file $tfile2 -order vtxyz mri_delete_chunk -dataset $tfile2 -chunk missing set rawy_args = ( $rawy_args $tfile2 ) *************** *** 308,325 **** mri_remap -file rawy -order tz -extents ${tdim}:${fakezdim} mri_remap -file scale -order tz -extents ${tdim}:${fakezdim} ! # Calculate the GLM ! echo "######## Calculating the GLM- this may take a while ########" echo "The command being used is:" ! echo mri_glm -v -ssqr -i rawy -p glm_results:images -scale scale \ ! $factoropts ! mri_glm -v -ssqr -i rawy -p glm_results:images -scale scale \ ! $factoropts ! echo '######## GLM complete; constructing results ########' # Map the GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i glm_results -f images.extent.v` ! mri_remap -file glm_results -order vxyz \ -extents ${vdim}:${xdim}:${ydim}:${zdim} # Grab interesting parts of the GLM. The output of the model consists --- 314,345 ---- mri_remap -file rawy -order tz -extents ${tdim}:${fakezdim} mri_remap -file scale -order tz -extents ${tdim}:${fakezdim} ! # Calculate the GLM over the raw (effect, group, and cross) factors ! echo "######## Calculating the first GLM- this may take a while ########" echo "The command being used is:" ! echo mri_glm -v -ssqr -i rawy -p raw_glm_results:images -o raw_residuals \ ! -scale scale $rawfactoropts ! mri_glm -v -ssqr -i rawy -p raw_glm_results:images -o raw_residuals \ ! -scale scale $rawfactoropts ! echo '######## First GLM complete ########' ! ! # Map this GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i raw_glm_results -f images.extent.v` ! mri_remap -file raw_glm_results -order vxyz \ ! -extents ${vdim}:${xdim}:${ydim}:${zdim} ! ! # Calculate the GLM over the subject factors ! echo "######## Calculating the second GLM- this may take a while ########" ! echo "The command being used is:" ! echo mri_glm -v -ssqr -i raw_residuals -p subj_glm_results:images \ ! -scale scale $subjfactoropts ! mri_glm -v -ssqr -i raw_residuals -p subj_glm_results:images \ ! -scale scale $subjfactoropts ! echo '######## Second GLM complete ########' # Map the GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i subj_glm_results -f images.extent.v` ! mri_remap -file subj_glm_results -order vxyz \ -extents ${vdim}:${xdim}:${ydim}:${zdim} # Grab interesting parts of the GLM. The output of the model consists *************** *** 329,366 **** # (nfactors+2) are the sums of squares- SSTO at offset (nfactors+2), # SSE at (nfactors+3), factor 0 SSR at (nfactors+4), and factor3 SSR # at (nfactors+5). ! mri_subset -d v -l 1 -s $groupfactor glm_results group_b ! mri_subset -d v -l 1 -s $crossfactor glm_results cross_b ! mri_subset -d v -l 1 -s $effectfactor glm_results effect_b ! @ ssto_offset = $nfactors + 2 ! mri_subset -d v -l 1 -s $ssto_offset glm_results ssto ! @ sse_offset = $nfactors + 3 ! mri_subset -d v -l 1 -s $sse_offset glm_results sse ! @ essqr_offset = $effectfactor + $nfactors + 3 ! mri_subset -d v -l 1 -s $essqr_offset glm_results ssb ! @ cssqr_offset = $crossfactor + $nfactors + 3 ! mri_subset -d v -l 1 -s $cssqr_offset glm_results ssc ! @ gssqr_offset = $groupfactor + $nfactors + 3 ! mri_subset -d v -l 1 -s $gssqr_offset glm_results ssg ! mri_subset -d v -l 1 -s 0 glm_results $homedir/mean echo "######## Generating final statistics ########" # The number of degrees of freedom in the sse is equal to the number of # input cases minus the number of fit parameters, or: ! # (nsubjs*nconds) - (nfactors + 1). ! # We use this to calculate MSE. ! @ sse_dof = ( $nsubjs * $nconds) - ( $nfactors + 1 ) mri_rpn_math -o mse '$1,'$sse_dof',/' sse # sssg is summed squared error for subject within group; we calculate it # from SSTO = SSR + SSE (with a check for negative values due to floating # point error). Degrees of freedom is the number of such columns in the # factor matrix, or: # $nsubjs - 2 . @ sssg_dof = $nsubjs - 2 ! mri_rpn_math -o sssg '$1,$2,$3,+,$4,+,$5,+,-,dup,0.0,swap,0.0,>,if_keep' \ ! ssto sse ssb ssc ssg mri_rpn_math -o mssg '$1,'$sssg_dof',/' sssg # Hop back home --- 349,394 ---- # (nfactors+2) are the sums of squares- SSTO at offset (nfactors+2), # SSE at (nfactors+3), factor 0 SSR at (nfactors+4), and factor3 SSR # at (nfactors+5). ! mri_subset -d v -l 1 -s $effectfactor raw_glm_results effect_b ! mri_subset -d v -l 1 -s $groupfactor raw_glm_results group_b ! mri_subset -d v -l 1 -s $crossfactor raw_glm_results cross_b ! @ ssto_offset = $nrawfactors + 2 ! mri_subset -d v -l 1 -s $ssto_offset raw_glm_results ssto ! @ sse_offset = $nrawfactors + 3 ! mri_subset -d v -l 1 -s $sse_offset raw_glm_results sse ! @ essqr_offset = $effectfactor + $nrawfactors + 3 ! mri_subset -d v -l 1 -s $essqr_offset raw_glm_results ssb ! @ gssqr_offset = $groupfactor + $nrawfactors + 3 ! mri_subset -d v -l 1 -s $gssqr_offset raw_glm_results ssg ! @ cssqr_offset = $crossfactor + $nrawfactors + 3 ! mri_subset -d v -l 1 -s $cssqr_offset raw_glm_results ssc ! @ subj_ssto_offset = $nsubjfactors + 2 ! mri_subset -d v -l 1 -s $subj_ssto_offset subj_glm_results subj_ssto ! @ subj_sse_offset = $nsubjfactors + 3 ! mri_subset -d v -l 1 -s $subj_sse_offset subj_glm_results subj_sse ! mri_subset -d v -l 1 -s 0 raw_glm_results $homedir/mean echo "######## Generating final statistics ########" # The number of degrees of freedom in the sse is equal to the number of # input cases minus the number of fit parameters, or: ! # (nsubjs*nconds) - (nrawfactors + 1). ! # We use this to calculate MSE (fixed effect model). ! @ sse_dof = ( $nsubjs * $nconds) - ( $nrawfactors + 1 ) mri_rpn_math -o mse '$1,'$sse_dof',/' sse + # We'll also need the degrees of freedom in the final sse (random + # effects model) after the second regression. + @ subj_sse_dof = ( $nsubjs * $nconds) - ( $nfactors + 1 ) + # sssg is summed squared error for subject within group; we calculate it # from SSTO = SSR + SSE (with a check for negative values due to floating # point error). Degrees of freedom is the number of such columns in the # factor matrix, or: # $nsubjs - 2 . @ sssg_dof = $nsubjs - 2 ! mri_rpn_math -o sssg '$1,$2,-,dup,0.0,swap,0.0,>,if_keep' \ ! subj_ssto subj_sse mri_rpn_math -o mssg '$1,'$sssg_dof',/' sssg # Hop back home *************** *** 371,385 **** # Following NKNW ch. 24, the relevant test statistic is MSB/MSAB where # factor B is the experimental condition (fixed) and factor A is the # subject identifier (random). This gives a pretty simple test statistic. mri_rpn_math -o $homedir/effect_f '$1,$2,/' $tmpdir/ssb $tmpdir/mssg ! mri_rpn_math -o $homedir/effect_p '$1,1,'$sse_dof',cf' $homedir/effect_f ! ! #mri_rpn_math -o $homedir/group_f '$1,$2,/' $tmpdir/ssg $tmpdir/mse ! #mri_rpn_math -o $homedir/group_p '$1,1,'$sssg_dof',cf' $homedir/group_f mri_rpn_math -o $homedir/cross_f '$1,$2,/' $tmpdir/ssc $tmpdir/mssg ! mri_rpn_math -o $homedir/cross_p '$1,1,'$sse_dof','$sssg_dof',+,cf' \ ! $homedir/cross_f # Clean up echo '######## Cleaning up #######' --- 399,417 ---- # Following NKNW ch. 24, the relevant test statistic is MSB/MSAB where # factor B is the experimental condition (fixed) and factor A is the # subject identifier (random). This gives a pretty simple test statistic. + # We do an additional small song-and-dance to cause the Pmaps to be saved + # in double precision. mri_rpn_math -o $homedir/effect_f '$1,$2,/' $tmpdir/ssb $tmpdir/mssg ! mri_type_convert -double $homedir/effect_f $tmpdir/tmp_doubles ! mri_rpn_math -o $homedir/effect_p '$1,1,'$sssg_dof',cf' $tmpdir/tmp_doubles mri_rpn_math -o $homedir/cross_f '$1,$2,/' $tmpdir/ssc $tmpdir/mssg ! mri_type_convert -double $homedir/cross_f $tmpdir/tmp_doubles ! mri_rpn_math -o $homedir/cross_p '$1,1,'$sssg_dof',cf' $tmpdir/tmp_doubles ! ! # Let's keep totalcounts as well; it's useful for seeing where the ! # samples were actually taken. ! mri_copy_dataset $tmpdir/totalcounts $homedir/totalcounts # Clean up echo '######## Cleaning up #######' Index: ./src/csh/Makefile *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/Makefile Wed Feb 13 19:01:16 2002 --- ./src/csh/Makefile Thu Apr 4 13:49:48 2002 *************** *** 19,25 **** MAKEFILES= Makefile SCRIPTFILES = FIASCO *.csh *.py SCRIPTSWITHDOCS = FIASCO groupcompare.csh merge_fisher.csh \ ! merge_detrend.csh false_discovery.py testingroup.csh include ../Makefile.pkg --- 19,26 ---- MAKEFILES= Makefile SCRIPTFILES = FIASCO *.csh *.py SCRIPTSWITHDOCS = FIASCO groupcompare.csh merge_fisher.csh \ ! merge_detrend.csh false_discovery.py testingroup.csh \ ! pooled_mean.csh pooled_stdv.csh include ../Makefile.pkg *************** *** 40,44 **** releaseprep: if test ! \( -f fiasco.local_orig.csh \); \ then cp fiasco.local.csh fiasco.local_orig.csh; \ ! else echo 'fiasco.local_orig.csh already exists' ; fi; chmod +xr $(SCRIPTFILES) --- 41,46 ---- releaseprep: if test ! \( -f fiasco.local_orig.csh \); \ then cp fiasco.local.csh fiasco.local_orig.csh; \ ! else echo 'fiasco.local_orig.csh already exists' ; fi chmod +xr $(SCRIPTFILES) + Index: ./src/csh/partialk.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/partialk.csh Wed Feb 13 19:01:18 2002 --- ./src/csh/partialk.csh Tue Mar 26 16:25:26 2002 *************** *** 27,37 **** # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '$Id: partialk.csh,v 1.5 2001/12/21 01:43:41 welling Exp $' # echo '#'`date`$0 if(! -d par) mkdir par ! @ newdim = 2 * $F_CLIP1_YCENTER partialk -reverse $F_PARTIALK_REVERSE -dim $newdim \ -input $1.mri -output $2.mri --- 27,41 ---- # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '$Id: partialk.csh,v 1.6 2002/03/26 21:25:33 welling Exp $' # echo '#'`date`$0 if(! -d par) mkdir par ! if ( $F_CLIP1_YCENTER != "center" ) then ! @ newdim = 2 * $F_CLIP1_YCENTER ! else ! @ newdim = `mri_printfield -input $1 -field "images.extent.y" -nofail` ! endif partialk -reverse $F_PARTIALK_REVERSE -dim $newdim \ -input $1.mri -output $2.mri Index: ./src/csh/pmaps.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/pmaps.csh Wed Feb 13 19:01:18 2002 --- ./src/csh/pmaps.csh Mon Apr 1 18:01:29 2002 *************** *** 1,5 **** #!/bin/csh -fx ! # fmaps.csh #/************************************************************ # * * # * Permission is hereby granted to any individual or * --- 1,5 ---- #!/bin/csh -fx ! # pmaps.csh #/************************************************************ # * * # * Permission is hereby granted to any individual or * *************** *** 29,35 **** # echo '#'`date`$0 ! echo '#$Id: pmaps.csh,v 1.1 2001/11/22 00:11:46 welling Exp $' # What shall we use for input? pushd $F_STAT_TPATH --- 29,35 ---- # echo '#'`date`$0 ! echo '#$Id: pmaps.csh,v 1.2 2002/04/01 23:01:41 welling Exp $' # What shall we use for input? pushd $F_STAT_TPATH *************** *** 62,76 **** mri_remap -file tmp_counts2 -order xyz mri_interp -d x -e $xdim -c tmp_counts2 tmp2_counts2 mri_interp -d y -e $ydim -c tmp2_counts2 tmp3_counts2 ! mri_rpn_math -o $p '$1,$2,$3,+,ct' $j tmp3_counts1 tmp3_counts2 mri_rpn_math -o ${p}_bar '$1,-1,*,$2,$3,+,ct' \ ! $j tmp3_counts1 tmp3_counts2 mri_destroy_dataset tmp_counts1 mri_destroy_dataset tmp2_counts1 mri_destroy_dataset tmp3_counts1 mri_destroy_dataset tmp_counts2 mri_destroy_dataset tmp2_counts2 mri_destroy_dataset tmp3_counts2 endif end --- 62,79 ---- mri_remap -file tmp_counts2 -order xyz mri_interp -d x -e $xdim -c tmp_counts2 tmp2_counts2 mri_interp -d y -e $ydim -c tmp2_counts2 tmp3_counts2 ! mri_type_convert -double $j tmp_dbl_tmap ! mri_rpn_math -o $p '$1,$2,$3,+,ct' \ ! tmp_dbl_tmap tmp3_counts1 tmp3_counts2 mri_rpn_math -o ${p}_bar '$1,-1,*,$2,$3,+,ct' \ ! tmp_dbl_tmap tmp3_counts1 tmp3_counts2 mri_destroy_dataset tmp_counts1 mri_destroy_dataset tmp2_counts1 mri_destroy_dataset tmp3_counts1 mri_destroy_dataset tmp_counts2 mri_destroy_dataset tmp2_counts2 mri_destroy_dataset tmp3_counts2 + mri_destroy_dataset tmp_dbl_tmap endif end Index: ./src/csh/pooled_mean.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/pooled_mean.csh Thu Apr 4 14:02:04 2002 --- ./src/csh/pooled_mean.csh Wed Mar 20 15:23:32 2002 *************** *** 0 **** --- 1,105 ---- + #! /usr/bin/csh -exf + # + # This script takes a series of Mean files as input (each of which + # also contains a counts chunk) and produces a single pooled Mean + # file as output. + # + + echo '#'`date` $0 + echo '# $Id: pooled_mean.csh,v 1.1 2002/03/20 20:24:28 welling Exp $' + + # Check for help command + if ( $#argv >= 1 ) then + if ( junk$argv[1] == "junk-help" ) then + if ( $#argv >= 2 ) then + scripthelp $0 $argv[2] + else + scripthelp $0 + endif + exit + endif + endif + if ( $#argv == 0 ) then + scripthelp $0 + exit + endif + + # Parse arguments + set outfile = "PooledMean" + set badarg = 0 + while ( junk$argv[1] =~ junk-* ) + switch ( $argv[1] ) + case '-out' : + set outfile = $argv[2]; shift argv; shift argv; breaksw; + default : + set badarg = 1; shift argv; breaksw; + endsw + if ( $#argv == 0 ) break; + end + if ( $badarg != 0 ) then + echo $0 ': invalid command line argument' + scripthelp $0 + exit -1 + endif + if ( $#argv == 0 ) then + echo $0 ': no mean files to pool!' + scripthelp $0 + exit -1 + endif + # + # Make some scratch space + # + if ( ${?F_TEMP} ) then + set tmpdir = ${F_TEMP}/pooled_mean_$$ + else + set tmpdir = ./pooled_mean_$$ + endif + if (! -e $tmpdir) mkdir $tmpdir + echo "Scratch directory will be " $tmpdir + + # Collect all means into a single file, indexed over t dimension + # We take some steps to restore the original dimensionality of + # the counts chunk. + set original_dims = `mri_printfield -input $1 -field counts.dimensions` + foreach fname ( $argv ) + mri_remap -file $fname -chunk counts -order zt + end + set fnamelist = ( $argv ) + mri_copy_dataset ${fnamelist[1]} $tmpdir/allmeans + shift fnamelist + while ( $#fnamelist ) + if ( $#fnamelist > 20 ) then + set these_args = ( $fnamelist[1-20] ) + set fnamelist = ( $fnamelist[21-] ) + else + set these_args = ( $fnamelist ) + set fnamelist = ( ) + endif + mri_paste -d t -out $tmpdir/tmp_paste $tmpdir/allmeans $these_args + mv $tmpdir/tmp_paste.mri $tmpdir/allmeans.mri + mv $tmpdir/tmp_paste.dat $tmpdir/allmeans.dat + end + foreach fname ( $argv ) + mri_remap -file $fname -chunk counts -order $original_dims + end + + # Build the pooled mean dataset + mri_copy_chunk -input $tmpdir/allmeans -output $tmpdir/allcounts \ + -chunk counts -to images -replace + mri_permute -input $tmpdir/allmeans -headerout $tmpdir/allmeans_p -order ztxy + mri_rpn_math -o $tmpdir/scaledmeans '$1,$2,*' \ + $tmpdir/allmeans_p $tmpdir/allcounts + mri_subsample -d t -e 1 -sum $tmpdir/scaledmeans $tmpdir/totalmeans + mri_subsample -d t -e 1 -sum $tmpdir/allcounts $tmpdir/totalcounts + mri_rpn_math -o $tmpdir/means_p '$1,$2,/' \ + $tmpdir/totalmeans $tmpdir/totalcounts + mri_permute -input $tmpdir/means_p -headerout $outfile -order xyzt + mri_copy_chunk -input $tmpdir/totalcounts -output $outfile -chunk images \ + -to counts -replace + mri_remap -file $outfile -chunk counts -order $original_dims + + # Clean up + rm -r $tmpdir + + + Index: ./src/csh/pooled_stdv.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/pooled_stdv.csh Thu Apr 4 14:02:10 2002 --- ./src/csh/pooled_stdv.csh Wed Mar 20 15:23:32 2002 *************** *** 0 **** --- 1,112 ---- + #! /usr/bin/csh -exf + # + # This script takes a series of Stdv files as input (each of which + # also contains a counts chunk) and produces a single pooled Stdv + # file as output. + # + echo '#'`date` $0 + echo '# $Id: pooled_stdv.csh,v 1.1 2002/03/20 20:24:28 welling Exp $' + + # Check for help command + if ( $#argv >= 1 ) then + if ( junk$argv[1] == "junk-help" ) then + if ( $#argv >= 2 ) then + scripthelp $0 $argv[2] + else + scripthelp $0 + endif + exit + endif + endif + if ( $#argv == 0 ) then + scripthelp $0 + exit + endif + + # Parse arguments + set outfile = "PooledStdv" + set badarg = 0 + while ( junk$argv[1] =~ junk-* ) + switch ( $argv[1] ) + case '-out' : + set outfile = $argv[2]; shift argv; shift argv; breaksw; + default : + set badarg = 1; shift argv; breaksw; + endsw + if ( $#argv == 0 ) break; + end + if ( $badarg != 0 ) then + echo $0 ': invalid command line argument' + scripthelp $0 + exit -1 + endif + if ( $#argv == 0 ) then + echo $0 ': no stdv files to pool!' + scripthelp $0 + exit -1 + endif + # + # Make some scratch space + # + if ( ${?F_TEMP} ) then + set tmpdir = ${F_TEMP}/pooled_stdv_$$ + else + set tmpdir = ./pooled_stdv_$$ + endif + if (! -e $tmpdir) mkdir $tmpdir + echo "Scratch directory will be " $tmpdir + + # Collect all stdvs into a single file, indexed over t dimension + # We take some steps to restore the original dimensionality of + # the counts chunk. + set original_dims = `mri_printfield -input $1 -field counts.dimensions` + foreach fname ( $argv ) + mri_remap -file $fname -chunk counts -order zt + end + set fnamelist = ( $argv ) + mri_copy_dataset ${fnamelist[1]} $tmpdir/allstdvs + shift fnamelist + while ( $#fnamelist ) + if ( $#fnamelist > 20 ) then + set these_args = ( $fnamelist[1-20] ) + set fnamelist = ( $fnamelist[21-] ) + else + set these_args = ( $fnamelist ) + set fnamelist = ( ) + endif + mri_paste -d t -out $tmpdir/tmp_paste $tmpdir/allstdvs $these_args + mv $tmpdir/tmp_paste.mri $tmpdir/allstdvs.mri + mv $tmpdir/tmp_paste.dat $tmpdir/allstdvs.dat + end + foreach fname ( $argv ) + mri_remap -file $fname -chunk counts -order $original_dims + end + + # Build the pooled stdv dataset + # + # numerator_terms is (n-1)s^2 for each of the component stdv's + # + # numerator is the sum of those terms + # denominator is (sum of n's) - 2 + # result is the square root of the ratio of these two + # + mri_copy_chunk -input $tmpdir/allstdvs -output $tmpdir/allcounts \ + -chunk counts -to images -replace + mri_permute -input $tmpdir/allstdvs -headerout $tmpdir/allstdvs_p -order ztxy + mri_rpn_math -o $tmpdir/numerator_terms '$1,dup,*,$2,1,-,*' \ + $tmpdir/allstdvs_p $tmpdir/allcounts + mri_subsample -d t -e 1 -sum $tmpdir/numerator_terms $tmpdir/numerator + mri_subsample -d t -e 1 -sum $tmpdir/allcounts $tmpdir/totalcounts + mri_rpn_math -o $tmpdir/denominator '$1,2,-' $tmpdir/totalcounts + mri_rpn_math -o $tmpdir/stdvs_p '$1,$2,/,sqrt' \ + $tmpdir/numerator $tmpdir/denominator + mri_permute -input $tmpdir/stdvs_p -headerout $outfile -order xyzt + mri_copy_chunk -input $tmpdir/totalcounts -output $outfile -chunk images \ + -to counts -replace + mri_remap -file $outfile -chunk counts -order $original_dims + + # Clean up + rm -r $tmpdir + + + Index: ./src/csh/spiral.defaults.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/spiral.defaults.csh Wed Feb 13 19:01:18 2002 --- ./src/csh/spiral.defaults.csh Thu Mar 21 19:03:01 2002 *************** *** 27,33 **** # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '$Id: spiral.defaults.csh,v 1.24 2001/12/21 01:45:26 welling Exp $' if(! -d in)mkdir in # # --- 27,33 ---- # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '$Id: spiral.defaults.csh,v 1.25 2002/03/22 00:04:00 welling Exp $' if(! -d in)mkdir in # # *************** *** 93,98 **** --- 93,111 ---- # adjustment parameters setenv F_OUTLIER_PARMS outlier.par # + # phase_lock.csh + # + # phase_lock header file + setenv F_PHLOCK_OUTPUT phase_lock + # unsmoothed estimates + setenv F_PHLOCK_RWPARMS phase_lock_raw.par + # smoothed estimates + setenv F_PHLOCK_PARMS phase_lock.par + # smoother bandwidth + setenv F_PHLOCK_BAND 3.0 + # algorithm: "weighted" or "center" + setenv F_PHLOCK_ALG "weighted" + # # # spiral.meanc.csh # Index: ./src/csh/spiral.despike.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/spiral.despike.csh Wed Feb 13 19:01:18 2002 --- ./src/csh/spiral.despike.csh Wed Mar 20 18:57:44 2002 *************** *** 28,185 **** # ************************************************************/ # echo '#'`date` $0 ! echo '#$Id: spiral.despike.csh,v 1.10 2002/02/05 00:26:21 bakalj Exp $' ! ! # ! # spk_cutoff is the cutoff threshold (in terms of (mag-median)/(IQR)) ! # for the body of the data (all but image 0) ! # spk_t0_cutoff is the cutoff threshold for image 0, in terms of ! # (orig-smoothed)/(IQR) ! # spk_t0_band is the smoother bandwidth used for image t0. ! # ! set spk_cutoff = $F_DESPIKE_CUTOFF ! set spk_t0_band = $F_DESPIKE_T0_BND ! set spk_t0_cutoff = $F_DESPIKE_T0_CUT set tdim = `mri_printfield -input $1 -field samples.extent.t` @ normimages = $tdim - 1 - @ median_offset = $normimages / 2 - @ q3_offset = $normimages / 4 - @ q1_offset = ( 3 * $normimages ) / 4 # Make a working copy, separating t=0 from the rest of the times mri_copy_chunk -input $1 -output data/raw_samples_all \ -chunk samples -to images -replace ! mri_subset -d t -l 1 -s 0 data/raw_samples_all data/raw_samples_t0 ! mri_subset -d t -l $normimages -s 1 data/raw_samples_all data/raw_samples mri_destroy_dataset data/raw_samples_all ! ! # ! # Prepare median and quartiles for magnitude for t!=0. ! # We clip off the first image for special processing. ! # ! mri_complex_to_scalar -m data/raw_samples data/raw_samples_m ! mri_permute -memlimit 50000000 -input data/raw_samples_m \ ! -headerout data/raw_samples_m_prm -order tvpscz ! mri_remap -file data/raw_samples_m_prm -order vtpscz ! mri_sort data/raw_samples_m_prm data/raw_samples_m_sorted ! mri_subset -d t -l 1 -s $median_offset data/raw_samples_m_sorted data/m_median ! mri_subset -d t -l 1 -s $q1_offset data/raw_samples_m_sorted data/m_q1 ! mri_subset -d t -l 1 -s $q3_offset data/raw_samples_m_sorted data/m_q3 ! mri_destroy_dataset data/raw_samples_m_sorted ! mri_destroy_dataset data/raw_samples_m_prm ! ! # ! # Prepare median for phase for t!=0. ! # We clip off the first image for special processing. ! # ! mri_complex_to_scalar -u data/raw_samples data/raw_samples_p ! mri_permute -memlimit 50000000 -input data/raw_samples_p \ ! -headerout data/raw_samples_p_prm -order tvpscz ! mri_destroy_dataset data/raw_samples_p ! mri_remap -file data/raw_samples_p_prm -order vtpscz ! mri_sort data/raw_samples_p_prm data/raw_samples_p_sorted ! mri_subset -d t -l 1 -s $median_offset data/raw_samples_p_sorted data/p_median ! mri_destroy_dataset data/raw_samples_p_prm ! mri_destroy_dataset data/raw_samples_p_sorted ! ! # Reconstruct complex median for t != 0 ! mri_remap -file data/m_median -order vpsczt ! mri_remap -file data/m_q1 -order vpsczt ! mri_remap -file data/m_q3 -order vpsczt ! mri_remap -file data/p_median -order vpsczt ! mri_rpn_math -o data/r_median '$1,$2,cos,*' data/m_median data/p_median ! mri_rpn_math -o data/i_median '$1,$2,sin,*' data/m_median data/p_median ! mri_destroy_dataset data/p_median ! mri_paste -d v -out data/median data/r_median data/i_median ! mri_destroy_dataset data/r_median ! mri_destroy_dataset data/i_median ! ! # ! # A given sample is considered to be a spike if its "trigger" value ! # (computed below) is larger than $spk_cutoff (defined above). ! # ! mri_rpn_math -o data/trigger '$1,$2,-,$3,$4,-,/,abs' \ ! data/raw_samples_m data/m_median data/m_q3 data/m_q1 ! mri_destroy_dataset data/raw_samples_m ! mri_remap -file data/trigger -order vpsczt ! mri_destroy_dataset data/m_median ! mri_rpn_math -o data/is_a_spike '$1,'$spk_cutoff',<' data/trigger ! mri_destroy_dataset data/trigger ! ! # ! # By substituting for the spike points, we make filtered values for ! # all times other than t=0. ! # ! mri_permute -memlimit 50000000 -input data/is_a_spike -headerout data/is_a_spike_tmp -order pscztv ! mri_paste -d v -out data/is_a_spike_complex_tmp data/is_a_spike_tmp data/is_a_spike_tmp ! mri_destroy_dataset data/is_a_spike_tmp ! mri_permute -memlimit 50000000 -input data/is_a_spike_complex_tmp -headerout data/is_a_spike_complex -order vpsczt ! mri_destroy_dataset data/is_a_spike_complex_tmp ! mri_rpn_math -o data/samples_filtered \ ! '$1,$2,$3,if_keep' \ ! data/raw_samples data/median data/is_a_spike_complex mri_destroy_dataset data/raw_samples - mri_destroy_dataset data/median - mri_destroy_dataset data/is_a_spike_complex - - # - # Special processing for t=0. We guess at where the spikes are by - # looking for rapid changes in magnitude, taking only the positive - # changes in this case. (Otherwise we end up picking up samples - # adjacent to the spike because of the smoothing operation). - # - mri_complex_to_scalar -m data/raw_samples_t0 data/raw_samples_m_t0 - mri_complex_to_scalar -u data/raw_samples_t0 data/raw_samples_p_t0 - mri_smooth -d p -b $spk_t0_band data/raw_samples_m_t0 data/smooth_m_t0 - mri_rpn_math -o data/trigger_t0 '$1,$2,-,$3,$4,-,/' \ - data/raw_samples_m_t0 data/smooth_m_t0 data/m_q3 data/m_q1 - mri_destroy_dataset data/m_q1 - mri_destroy_dataset data/m_q3 - mri_destroy_dataset data/smooth_m_t0 - mri_rpn_math -o data/is_a_spike_t0 '$1,'$spk_t0_cutoff',<' data/trigger_t0 - mri_paste -out data/is_a_spike_all -d t data/is_a_spike_t0 data/is_a_spike - mri_destroy_dataset data/is_a_spike - mri_destroy_dataset data/trigger_t0 - - # - # Through clever use of a triangular filter we happen to have handy, - # we'll produce a set of replacement values for the spikes which are - # the average of the two adjacent points. - # - mri_smooth -d p -b 2 -kernel triangular \ - data/raw_samples_m_t0 data/triangle_m_t0 - mri_smooth -d p -b 2 -kernel triangular \ - data/raw_samples_p_t0 data/triangle_p_t0 - mri_rpn_math -o data/interp_m_t0 '$1,$2,0.5,*,-,2,*' \ - data/triangle_m_t0 data/raw_samples_m_t0 - mri_destroy_dataset data/triangle_m_t0 - mri_destroy_dataset data/raw_samples_m_t0 - mri_rpn_math -o data/interp_p_t0 '$1,$2,0.5,*,-,2,*' \ - data/triangle_p_t0 data/raw_samples_p_t0 - mri_destroy_dataset data/triangle_p_t0 - mri_destroy_dataset data/raw_samples_p_t0 ! # ! # Substitute these replacement values where appropriate ! # ! mri_rpn_math -o data/r_interp_t0 '$1,$2,cos,*' \ ! data/interp_m_t0 data/interp_p_t0 ! mri_rpn_math -o data/i_interp_t0 '$1,$2,sin,*' \ ! data/interp_m_t0 data/interp_p_t0 ! mri_destroy_dataset data/interp_m_t0 ! mri_destroy_dataset data/interp_p_t0 ! mri_paste -d v -out data/interp_t0 data/r_interp_t0 data/i_interp_t0 ! mri_destroy_dataset data/r_interp_t0 ! mri_destroy_dataset data/i_interp_t0 ! mri_interp -c -d v -e 2 data/is_a_spike_t0 data/is_a_spike_t0_complex ! mri_destroy_dataset data/is_a_spike_t0 ! mri_rpn_math -o data/samples_filtered_t0 \ ! '$1,$2,$3,if_keep' \ ! data/raw_samples_t0 data/interp_t0 data/is_a_spike_t0_complex mri_destroy_dataset data/raw_samples_t0 - mri_destroy_dataset data/interp_t0 - mri_destroy_dataset data/is_a_spike_t0_complex # Produce an output dataset mri_paste -out data/samples_filtered_all -d t \ --- 28,74 ---- # ************************************************************/ # echo '#'`date` $0 ! echo '#$Id: spiral.despike.csh,v 1.11 2002/03/20 23:59:12 welling Exp $' set tdim = `mri_printfield -input $1 -field samples.extent.t` @ normimages = $tdim - 1 # Make a working copy, separating t=0 from the rest of the times mri_copy_chunk -input $1 -output data/raw_samples_all \ -chunk samples -to images -replace ! mri_type_convert -float data/raw_samples_all data/raw_samples_all_float mri_destroy_dataset data/raw_samples_all ! mri_subset -d t -l 1 -s 0 data/raw_samples_all_float data/raw_samples_t0 ! mri_subset -d t -l $normimages -s 1 data/raw_samples_all_float data/raw_samples ! mri_destroy_dataset data/raw_samples_all_float ! ! # Despike the bulk of the data by comparing across times. ! # This produces some files that we want as side effects. ! despike_timewise.csh data/raw_samples data/samples_filtered data/spike_count ! ! if (! -f data/raw_samples_m_median.mri) then ! echo "despike_timewise.csh failed to produce magnitude median!" ! exit -1 ! endif ! if (! -f data/raw_samples_m_q1.mri) then ! echo "despike_timewise.csh failed to produce magnitude q1!" ! exit -1 ! endif ! if (! -f data/raw_samples_m_q3.mri) then ! echo "despike_timewise.csh failed to produce magnitude q3!" ! exit -1 ! endif ! if (! -f data/raw_samples_p_median.mri) then ! echo "despike_timewise.csh failed to produce phase median!" ! exit -1 ! endif mri_destroy_dataset data/raw_samples ! # Despike image 0 by comparing across samples. ! despike_sampwise.csh data/raw_samples_t0 \ ! data/raw_samples_m_q1 data/raw_samples_m_q3 \ ! data/samples_filtered_t0 data/spike_count_t0 mri_destroy_dataset data/raw_samples_t0 # Produce an output dataset mri_paste -out data/samples_filtered_all -d t \ *************** *** 193,207 **** # Produce summary information if(! -d par) mkdir par ! mri_subsample -sum -d p -e 1 data/is_a_spike_all data/spike_count set step = `depath.csh $0` echo '##Format: order:index_tz type:raw names:(spike_count)' \ > par/$F_DESPIKE_PARMS mri_rpn_math -o data/junk '0,$t,$z,$1,1,if_print_3' \ ! data/spike_count >> par/$F_DESPIKE_PARMS echo "$step par/${F_DESPIKE_PARMS}" >> $F_SUMM_INPUT mri_destroy_dataset data/junk ! mri_destroy_dataset data/spike_count - # Clean up - mri_destroy_dataset data/is_a_spike_all --- 82,96 ---- # Produce summary information if(! -d par) mkdir par ! mri_paste -out data/spike_count_all data/spike_count_t0 data/spike_count set step = `depath.csh $0` + mri_destroy_dataset data/spike_count + mri_destroy_dataset data/spike_count_t0 echo '##Format: order:index_tz type:raw names:(spike_count)' \ > par/$F_DESPIKE_PARMS mri_rpn_math -o data/junk '0,$t,$z,$1,1,if_print_3' \ ! data/spike_count_all >> par/$F_DESPIKE_PARMS echo "$step par/${F_DESPIKE_PARMS}" >> $F_SUMM_INPUT mri_destroy_dataset data/junk ! mri_destroy_dataset data/spike_count_all Index: ./src/csh/testingroup.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/testingroup.csh Wed Feb 13 19:01:19 2002 --- ./src/csh/testingroup.csh Thu Apr 4 13:48:07 2002 *************** *** 3,9 **** # This script does between-group comparisons. # echo '#'`date` $0 ! echo '# $Id: testingroup.csh,v 1.7 2002/01/10 18:34:06 bakalj Exp $' # Check for help command if ( $#argv >= 1 ) then --- 3,9 ---- # This script does between-group comparisons. # echo '#'`date` $0 ! echo '# $Id: testingroup.csh,v 1.8 2002/04/04 18:48:09 welling Exp $' # Check for help command if ( $#argv >= 1 ) then *************** *** 152,181 **** mri_rpn_math -o rawfactor1 '$c,2,*,1,-' conds_x_subjects set which_subj = 0 ! set factorlist = ( rawfactor1 ) set effectfactor = 1 foreach subj ($args) if ( $subj == $lastsubj ) break mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,1.0,'$nsubjs',/,-,$2,*' \ conds_x_subjects rawfactor1 ! set factorlist = ( $factorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end ! @ nfactors = ${#factorlist} ! echo "nfactors is " $nfactors ! echo "unscaled factors are " $factorlist echo "The factor of interest is number " $effectfactor # Remap the factors into the order needed for GLM ! foreach fname ( $factorlist ) mri_remap -file $fname -order t -extents ${tdim} end # Make the factor part of the operand string for mri_glm ! set factoropts = ( ) ! foreach fname ( $factorlist ) ! set factoropts = ( $factoropts -f ${fname}:images ) end # Build the scaling matrix and the vector of mean values. There are --- 152,190 ---- mri_rpn_math -o rawfactor1 '$c,2,*,1,-' conds_x_subjects set which_subj = 0 ! set rawfactorlist = ( rawfactor1 ) set effectfactor = 1 + set subjfactorlist = ( ) foreach subj ($args) if ( $subj == $lastsubj ) break mri_rpn_math -o normsubjfactor$which_subj \ '$s,'${which_subj}',==,1.0,'$nsubjs',/,-,$2,*' \ conds_x_subjects rawfactor1 ! set subjfactorlist = ( $subjfactorlist normsubjfactor${which_subj} ) @ which_subj = $which_subj + 1 end ! @ nrawfactors = ${#rawfactorlist} ! @ nsubjfactors = ${#subjfactorlist} ! @ nfactors = $nrawfactors + $nsubjfactors ! echo "number of fixed effect factors is " $nrawfactors ! echo "number of factors associated with individual subjects is " $nsubjfactors ! echo "total nfactors is " $nfactors ! echo "unscaled factors are " $rawfactorlist $subjfactorlist echo "The factor of interest is number " $effectfactor # Remap the factors into the order needed for GLM ! foreach fname ( $rawfactorlist $subjfactorlist ) mri_remap -file $fname -order t -extents ${tdim} end # Make the factor part of the operand string for mri_glm ! set rawfactoropts = ( ) ! foreach fname ( $rawfactorlist ) ! set rawfactoropts = ( $rawfactoropts -f ${fname}:images ) ! end ! set subjfactoropts = ( ) ! foreach fname ( $subjfactorlist ) ! set subjfactoropts = ( $subjfactoropts -f ${fname}:images ) end # Build the scaling matrix and the vector of mean values. There are *************** *** 270,308 **** mri_remap -file rawy -order tz -extents ${tdim}:${fakezdim} mri_remap -file scale -order tz -extents ${tdim}:${fakezdim} ! # Calculate the GLM ! echo "######## Calculating the GLM- this may take a while ########" echo "The command being used is:" ! echo mri_glm -v -ssqr -i rawy -p glm_results:images -scale scale \ ! $factoropts ! mri_glm -v -ssqr -i rawy -p glm_results:images -scale scale \ ! $factoropts ! echo '######## GLM complete; constructing results ########' # Map the GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i glm_results -f images.extent.v` ! mri_remap -file glm_results -order vxyz \ -extents ${vdim}:${xdim}:${ydim}:${zdim} # Grab interesting parts of the GLM. The output of the model consists # of (nfactors+1) estimated b's, followed by the associated summed squares. ! mri_subset -d v -l 1 -s $effectfactor glm_results effect_b ! @ ssto_offset = $nfactors + 2 ! mri_subset -d v -l 1 -s $ssto_offset glm_results ssto ! @ sse_offset = $nfactors + 3 ! mri_subset -d v -l 1 -s $sse_offset glm_results sse ! @ essqr_offset = $effectfactor + $nfactors + 3 ! mri_subset -d v -l 1 -s $essqr_offset glm_results ssb ! mri_subset -d v -l 1 -s 0 glm_results $homedir/mean echo "######## Generating final statistics ########" ! # The number of degrees of freedom in the sse is equal to the number of # input cases minus the number of fit parameters, or: # (nsubjs*nconds) - (nfactors + 1). ! # We use this to calculate MSE. ! @ sse_dof = ( $nsubjs * $nconds) - ( $nfactors + 1 ) ! mri_rpn_math -o mse '$1,'$sse_dof',/' sse # sssg is summed squared error for subject within group; we calculate it # from SSTO = SSR + SSE (with a check for negative values due to floating --- 279,335 ---- mri_remap -file rawy -order tz -extents ${tdim}:${fakezdim} mri_remap -file scale -order tz -extents ${tdim}:${fakezdim} ! # Calculate the GLM over the raw (effect) factors ! echo "######## Calculating the first GLM- this may take a while ########" echo "The command being used is:" ! echo mri_glm -v -ssqr -i rawy -p raw_glm_results:images -o raw_residuals \ ! -scale scale $rawfactoropts ! mri_glm -v -ssqr -i rawy -p raw_glm_results:images -o raw_residuals \ ! -scale scale $rawfactoropts ! echo '######## First GLM complete ########' # Map the GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i raw_glm_results -f images.extent.v` ! mri_remap -file raw_glm_results -order vxyz \ ! -extents ${vdim}:${xdim}:${ydim}:${zdim} ! ! # Calculate the GLM over the subject factors ! echo "######## Calculating the second GLM- this may take a while ########" ! echo "The command being used is:" ! echo mri_glm -v -ssqr -i raw_residuals -p subj_glm_results:images \ ! -scale scale $subjfactoropts ! mri_glm -v -ssqr -i raw_residuals -p subj_glm_results:images \ ! -scale scale $subjfactoropts ! echo '######## Second GLM complete; constructing results ########' ! ! # Map the GLM results back to appropriate dimensions ! set vdim = `mri_printfield -i subj_glm_results -f images.extent.v` ! mri_remap -file subj_glm_results -order vxyz \ -extents ${vdim}:${xdim}:${ydim}:${zdim} # Grab interesting parts of the GLM. The output of the model consists # of (nfactors+1) estimated b's, followed by the associated summed squares. ! mri_subset -d v -l 1 -s $effectfactor raw_glm_results effect_b ! @ ssto_offset = $nrawfactors + 2 ! mri_subset -d v -l 1 -s $ssto_offset raw_glm_results ssto ! @ sse_offset = $nrawfactors + 3 ! mri_subset -d v -l 1 -s $sse_offset raw_glm_results sse ! @ essqr_offset = $effectfactor + $nrawfactors + 3 ! mri_subset -d v -l 1 -s $essqr_offset raw_glm_results ssb ! @ subj_ssto_offset = $nsubjfactors + 2 ! mri_subset -d v -l 1 -s $subj_ssto_offset subj_glm_results subj_ssto ! @ subj_sse_offset = $nsubjfactors + 3 ! mri_subset -d v -l 1 -s $subj_sse_offset subj_glm_results subj_sse ! mri_subset -d v -l 1 -s 0 raw_glm_results $homedir/mean echo "######## Generating final statistics ########" ! # The number of degrees of freedom in the final sse is equal to the number of # input cases minus the number of fit parameters, or: # (nsubjs*nconds) - (nfactors + 1). ! # We use this to calculate the final MSE. ! @ subj_sse_dof = ( $nsubjs * $nconds) - ( $nfactors + 1 ) ! mri_rpn_math -o subj_mse '$1,'$subj_sse_dof',/' subj_sse # sssg is summed squared error for subject within group; we calculate it # from SSTO = SSR + SSE (with a check for negative values due to floating *************** *** 310,317 **** # factor matrix, or: # $nsubjs - 1 . @ sssg_dof = $nsubjs - 1 ! mri_rpn_math -o sssg '$1,$2,$3,+,-,dup,0.0,swap,0.0,>,if_keep' \ ! ssto sse ssb mri_rpn_math -o mssg '$1,'$sssg_dof',/' sssg # Hop back home --- 337,344 ---- # factor matrix, or: # $nsubjs - 1 . @ sssg_dof = $nsubjs - 1 ! mri_rpn_math -o sssg '$1,$2,-,dup,0.0,swap,0.0,>,if_keep' \ ! subj_ssto subj_sse mri_rpn_math -o mssg '$1,'$sssg_dof',/' sssg # Hop back home *************** *** 322,330 **** # Following NKNW ch. 24, the relevant test statistic is MSB/MSAB where # factor B is the experimental condition (fixed) and factor A is the # subject identifier (random). This gives a pretty simple test statistic. ! @ ndf = 1 mri_rpn_math -o $homedir/effect_f '$1,$2,/' $tmpdir/ssb $tmpdir/mssg ! mri_rpn_math -o $homedir/effect_p '$1,1,'$sse_dof',cf' $homedir/effect_f # Clean up echo '######## Cleaning up #######' --- 349,363 ---- # Following NKNW ch. 24, the relevant test statistic is MSB/MSAB where # factor B is the experimental condition (fixed) and factor A is the # subject identifier (random). This gives a pretty simple test statistic. ! # We do an additional small song-and-dance to cause the Pmaps to be saved ! # in double precision. mri_rpn_math -o $homedir/effect_f '$1,$2,/' $tmpdir/ssb $tmpdir/mssg ! mri_type_convert -double $homedir/effect_f $tmpdir/tmp_doubles ! mri_rpn_math -o $homedir/effect_p '$1,1,'$sssg_dof',cf' $tmpdir/tmp_doubles ! ! # Let's keep totalcounts as well; it's useful for seeing where the ! # samples were actually taken. ! mri_copy_dataset $tmpdir/totalcounts $homedir/totalcounts # Clean up echo '######## Cleaning up #######' Index: ./src/csh/ts.defaults.csh *** /TOIL-U4/fiasco/Fiasco5.1/src/csh/ts.defaults.csh Wed Feb 13 19:01:19 2002 --- ./src/csh/ts.defaults.csh Thu Mar 21 19:02:48 2002 *************** *** 27,33 **** # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '#$Id: ts.defaults.csh,v 1.30 2001/12/21 01:45:26 welling Exp $' if(! -d in)mkdir in # # --- 27,33 ---- # * Original programming by Bill Eddy * # ************************************************************/ # ! echo '#$Id: ts.defaults.csh,v 1.32 2002/03/22 00:03:46 welling Exp $' if(! -d in)mkdir in # # *************** *** 98,103 **** --- 98,110 ---- setenv F_READER_OPTS "" # # + # epi.despike.csh + # + setenv F_DESPIKE_OUTPUT despike + setenv F_DESPIKE_PARMS spike_count.par + setenv F_DESPIKE_CUTOFF 4.0 + # + # # ts.baseline.csh # # baseline header file *************** *** 172,177 **** --- 179,197 ---- # deghost filter bandwidth setenv F_DEGHOST2_BAND 8.0 # + # phase_lock.csh + # + # phase_lock header file + setenv F_PHLOCK_OUTPUT phase_lock + # unsmoothed estimates + setenv F_PHLOCK_RWPARMS phase_lock_raw.par + # smoothed estimates + setenv F_PHLOCK_PARMS phase_lock.par + # smoother bandwidth + setenv F_PHLOCK_BAND 3.0 + # algorithm: "weighted" or "center" + setenv F_PHLOCK_ALG "weighted" + # # ts.meanc.csh # # meanc header file Index: ./src/fmri/fmri.c *** /TOIL-U4/fiasco/Fiasco5.1/src/fmri/fmri.c Wed Feb 13 19:01:11 2002 --- ./src/fmri/fmri.c Thu Mar 28 18:15:43 2002 *************** *** 54,60 **** #include "fmri.h" #include "stdcrg.h" ! static char rcsid[] = "$Id: fmri.c,v 1.9 2002/01/29 21:26:00 welling Exp $"; /* Get or create a vector of indicators for missing images */ unsigned char **get_missing( MRI_Dataset* DS ) --- 54,60 ---- #include "fmri.h" #include "stdcrg.h" ! static char rcsid[] = "$Id: fmri.c,v 1.10 2002/03/28 23:15:28 welling Exp $"; /* Get or create a vector of indicators for missing images */ unsigned char **get_missing( MRI_Dataset* DS ) *************** *** 117,124 **** Abort("Input dataset missing chunk has no datatype!\n"); if (strcmp( mri_get_string( DS, "missing.datatype" ), "uint8" )) Abort("Input dataset missing chunk has wrong datatype!\n"); ! missing_dz = mri_get_int( DS, "images.extent.z" ); ! missing_dt = mri_get_int( DS, "images.extent.t" ); } /* If both are present, are they consistent? */ --- 117,124 ---- Abort("Input dataset missing chunk has no datatype!\n"); if (strcmp( mri_get_string( DS, "missing.datatype" ), "uint8" )) Abort("Input dataset missing chunk has wrong datatype!\n"); ! missing_dz = mri_get_int( DS, "missing.extent.z" ); ! missing_dt = mri_get_int( DS, "missing.extent.t" ); } /* If both are present, are they consistent? */ Index: ./src/fmri/fshrot3d.c *** /TOIL-U4/fiasco/Fiasco5.1/src/fmri/fshrot3d.c Wed Feb 13 19:01:11 2002 --- ./src/fmri/fshrot3d.c Thu Apr 4 13:51:00 2002 *************** *** 55,61 **** * -Note that the routine expects data presented in z-fastest order! */ ! static char rcsid[] = "$Id: fshrot3d.c,v 1.15 2001/01/03 21:52:48 welling Exp $"; /* How bad does cancellation have to get before we consider a quat bad? */ #define CANCELLATION_TOL 0.01 --- 55,61 ---- * -Note that the routine expects data presented in z-fastest order! */ ! static char rcsid[] = "$Id: fshrot3d.c,v 1.16 2002/04/04 18:50:33 welling Exp $"; /* How bad does cancellation have to get before we consider a quat bad? */ #define CANCELLATION_TOL 0.01 *************** *** 977,987 **** if (debug) fprintf(stderr,"shear_x %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx % 2) nx_mod= nx-1; else nx_mod= nx; ! if (ny % 2) ny_mod= ny-1; else ny_mod= ny; ! if (nz % 2) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ --- 977,987 ---- if (debug) fprintf(stderr,"shear_x %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx>1 && (nx % 2)) nx_mod= nx-1; else nx_mod= nx; ! if (ny>1 && (ny % 2)) ny_mod= ny-1; else ny_mod= ny; ! if (nz>1 && (nz % 2)) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ *************** *** 1055,1065 **** if (debug) fprintf(stderr,"shear_y %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx % 2) nx_mod= nx-1; else nx_mod= nx; ! if (ny % 2) ny_mod= ny-1; else ny_mod= ny; ! if (nz % 2) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ --- 1055,1065 ---- if (debug) fprintf(stderr,"shear_y %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx>1 && (nx % 2)) nx_mod= nx-1; else nx_mod= nx; ! if (ny>1 && (ny % 2)) ny_mod= ny-1; else ny_mod= ny; ! if (nz>1 && (nz % 2)) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ *************** *** 1132,1142 **** if (debug) fprintf(stderr,"shear_z %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx % 2) nx_mod= nx-1; else nx_mod= nx; ! if (ny % 2) ny_mod= ny-1; else ny_mod= ny; ! if (nz % 2) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ --- 1132,1142 ---- if (debug) fprintf(stderr,"shear_z %f %f %f\n",a,b,delta); /* Adjustments for odd data dimensions */ ! if (nx>1 && (nx % 2)) nx_mod= nx-1; else nx_mod= nx; ! if (ny>1 && (ny % 2)) ny_mod= ny-1; else ny_mod= ny; ! if (nz>1 && (nz % 2)) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ *************** *** 1608,1618 **** dx, dy, dz); /* Adjustments for odd data dimensions */ ! if (nx % 2) nx_mod= nx-1; else nx_mod= nx; ! if (ny % 2) ny_mod= ny-1; else ny_mod= ny; ! if (nz % 2) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ --- 1608,1618 ---- dx, dy, dz); /* Adjustments for odd data dimensions */ ! if (nx>1 && (nx % 2)) nx_mod= nx-1; else nx_mod= nx; ! if (ny>1 && (ny % 2)) ny_mod= ny-1; else ny_mod= ny; ! if (nz>1 && (nz % 2)) nz_mod= nz-1; else nz_mod= nz; /* Step counter */ Index: ./src/fmri/glm.c *** /TOIL-U4/fiasco/Fiasco5.1/src/fmri/glm.c Wed Feb 13 19:01:11 2002 --- ./src/fmri/glm.c Tue Apr 2 19:47:16 2002 *************** *** 79,99 **** static char err_buf[512]; /* not used for all messages */ static int mach_prec_initialized= 0; ! static float mach_precision; static int old_nfactors= 0; static int old_nobs= 0; static int old_complex_flag= 0; static int last_fit_valid= 0; /* Various scratch areas */ ! static float* factors_ftn= NULL; ! static float* u_ftn= NULL; ! static float* v_ftn= NULL; ! static float* singular= NULL; ! static float* inv_singular= NULL; ! static float* residuals= NULL; ! static float* covariance_ftn= NULL; ! static float* rwork= NULL; static int lwork= 0; static int gbl_complex_flag= 0; --- 79,99 ---- static char err_buf[512]; /* not used for all messages */ static int mach_prec_initialized= 0; ! static double mach_precision; static int old_nfactors= 0; static int old_nobs= 0; static int old_complex_flag= 0; static int last_fit_valid= 0; /* Various scratch areas */ ! static double* factors_ftn= NULL; ! static double* u_ftn= NULL; ! static double* v_ftn= NULL; ! static double* singular= NULL; ! static double* inv_singular= NULL; ! static double* residuals= NULL; ! static double* covariance_ftn= NULL; ! static double* rwork= NULL; static int lwork= 0; static int gbl_complex_flag= 0; *************** *** 158,164 **** return err_txt; } ! static void flip_to_ftn( float* in, float* out, int xdim, int ydim ) { int i; int j; --- 158,164 ---- return err_txt; } ! static void flip_to_ftn( double* in, double* out, int xdim, int ydim ) { int i; int j; *************** *** 167,173 **** for (j=0; jdlim) dlim= data_in[i]; ! #ifdef __GNUC__ ! dlim *= (float)sqrt((float)nsamples) * mach_precision; ! #else ! dlim *= sqrtf((float)nsamples) * mach_precision; ! #endif for (i=0; idlim) dlim= data_in[i]; ! dlim *= sqrt((double)nsamples) * mach_precision; for (i=0; i MAXFLOAT) ! { ! float_buf[i] = MAXFLOAT; ! conv_error = TRUE; ! } ! else float_buf[i] = (float) double_buf[i]; mri_set_chunk(ds, key, size, offset, MRI_FLOAT, float_buf); break; case MRI_FLOAT64: --- 1888,1910 ---- case MRI_FLOAT32: float_buf = (float *) GetBuffer(ds, size*sizeof(float)); for (i = 0; i < size; ++i) ! if (finite(double_buf[i])) { ! if (double_buf[i] < -MAXFLOAT) ! { ! float_buf[i] = -MAXFLOAT; ! conv_error = TRUE; ! } ! else if (double_buf[i] > MAXFLOAT) ! { ! float_buf[i] = MAXFLOAT; ! conv_error = TRUE; ! } ! else float_buf[i] = (float) double_buf[i]; ! } ! else { ! /* We'll let the compiler translate the NaN or Inf */ ! float_buf[i]= (float)double_buf[i]; ! } mri_set_chunk(ds, key, size, offset, MRI_FLOAT, float_buf); break; case MRI_FLOAT64: Index: ./src/mri_anova/mri_anova.c *** /TOIL-U4/fiasco/Fiasco5.1/src/mri_anova/mri_anova.c Wed Feb 13 19:01:23 2002 --- ./src/mri_anova/mri_anova.c Wed Mar 27 19:58:06 2002 *************** *** 72,78 **** /* Notes- */ ! static char rcsid[] = "$Id: mri_anova.c,v 1.11 2000/11/02 20:02:15 welling Exp $"; #define MAX_SOURCES 10 #define CHUNKSIZE 4096 --- 72,78 ---- /* Notes- */ ! static char rcsid[] = "$Id: mri_anova.c,v 1.12 2002/03/28 00:57:57 welling Exp $"; #define MAX_SOURCES 10 #define CHUNKSIZE 4096 *************** *** 88,94 **** typedef struct out_file_struct { char* source; MRI_Dataset* dset; ! float* chunk; int chunk_offset; int file_offset; } OutFileStruct; --- 88,94 ---- typedef struct out_file_struct { char* source; MRI_Dataset* dset; ! double* chunk; int chunk_offset; int file_offset; } OutFileStruct; *************** *** 106,116 **** static sp_SplitRec* split_table= NULL; static OutFileStruct* out_table= NULL; static OutFileStruct* mse_out= NULL; ! static float* factors_unpacked= NULL; ! static float* factors= NULL; ! static float* factor_means= NULL; ! static float* tseries= NULL; ! static float* params= NULL; /* This routine generates an approprate dataset name from a * source name. --- 106,116 ---- static sp_SplitRec* split_table= NULL; static OutFileStruct* out_table= NULL; static OutFileStruct* mse_out= NULL; ! static double* factors_unpacked= NULL; ! static double* factors= NULL; ! static double* factor_means= NULL; ! static double* tseries= NULL; ! static double* params= NULL; /* This routine generates an approprate dataset name from a * source name. *************** *** 138,145 **** /* Build a prototype output ds in the first table element */ out= out_table; out->source= source_table[0].name; ! if (!(out->chunk= (float*)malloc(CHUNKSIZE*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",progname,CHUNKSIZE); out->chunk_offset= 0; out->file_offset= 0; out->dset= mri_open_dataset( fname_from_source(out->source), MRI_WRITE ); --- 138,145 ---- /* Build a prototype output ds in the first table element */ out= out_table; out->source= source_table[0].name; ! if (!(out->chunk= (double*)malloc(CHUNKSIZE*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",progname,CHUNKSIZE); out->chunk_offset= 0; out->file_offset= 0; out->dset= mri_open_dataset( fname_from_source(out->source), MRI_WRITE ); *************** *** 169,176 **** for (i=1; isource= source_table[i].name; ! if (!(out->chunk= (float*)malloc(CHUNKSIZE*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",progname,CHUNKSIZE); out->chunk_offset= 0; out->file_offset= 0; out->dset = mri_copy_dataset( fname_from_source(out->source), --- 169,176 ---- for (i=1; isource= source_table[i].name; ! if (!(out->chunk= (double*)malloc(CHUNKSIZE*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",progname,CHUNKSIZE); out->chunk_offset= 0; out->file_offset= 0; out->dset = mri_copy_dataset( fname_from_source(out->source), *************** *** 183,190 **** Abort("%s: unable to allocate %d bytes!\n", progname, sizeof(OutFileStruct)); mse_out->source= "MSE"; ! if (!(mse_out->chunk= (float*)malloc(CHUNKSIZE*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",progname,CHUNKSIZE); mse_out->chunk_offset= 0; mse_out->file_offset= 0; mse_out->dset= mri_copy_dataset( "MSE", out_table[0].dset ); --- 183,190 ---- Abort("%s: unable to allocate %d bytes!\n", progname, sizeof(OutFileStruct)); mse_out->source= "MSE"; ! if (!(mse_out->chunk= (double*)malloc(CHUNKSIZE*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",progname,CHUNKSIZE); mse_out->chunk_offset= 0; mse_out->file_offset= 0; mse_out->dset= mri_copy_dataset( "MSE", out_table[0].dset ); *************** *** 200,206 **** out= out_table+i; if (out->chunk_offset) { mri_set_chunk(out->dset, "images", out->chunk_offset, ! out->file_offset, MRI_FLOAT, out->chunk); out->file_offset += out->chunk_offset; out->chunk_offset= 0; } --- 200,206 ---- out= out_table+i; if (out->chunk_offset) { mri_set_chunk(out->dset, "images", out->chunk_offset, ! out->file_offset, MRI_DOUBLE, out->chunk); out->file_offset += out->chunk_offset; out->chunk_offset= 0; } *************** *** 208,214 **** if (mse_out) { if (mse_out->chunk_offset) { mri_set_chunk(mse_out->dset, "images", mse_out->chunk_offset, ! mse_out->file_offset, MRI_FLOAT, mse_out->chunk); mse_out->file_offset += mse_out->chunk_offset; mse_out->chunk_offset= 0; } --- 208,214 ---- if (mse_out) { if (mse_out->chunk_offset) { mri_set_chunk(mse_out->dset, "images", mse_out->chunk_offset, ! mse_out->file_offset, MRI_DOUBLE, mse_out->chunk); mse_out->file_offset += mse_out->chunk_offset; mse_out->chunk_offset= 0; } *************** *** 320,342 **** static void allocate_memory() { ! if (!(factors_unpacked= (float*)malloc(nimages*nsources*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",nimages*nsources); ! if (!(factors= (float*)malloc(nimages*nsources*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",nimages*nsources); ! if (!(factor_means= (float*)malloc(nsources*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",nimages*nsources); ! if (!(tseries= (float*)malloc(nimages*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",nimages); ! if (!(params= (float*)malloc(nparams*sizeof(float)))) ! Abort("%s: unable to allocate %d floats!\n",nparams); } ! static void build_factor_matrix(int z, float* fac) { int t; int s; --- 320,342 ---- static void allocate_memory() { ! if (!(factors_unpacked= (double*)malloc(nimages*nsources*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",nimages*nsources); ! if (!(factors= (double*)malloc(nimages*nsources*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",nimages*nsources); ! if (!(factor_means= (double*)malloc(nsources*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",nimages*nsources); ! if (!(tseries= (double*)malloc(nimages*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",nimages); ! if (!(params= (double*)malloc(nparams*sizeof(double)))) ! Abort("%s: unable to allocate %d doubles!\n",nparams); } ! static void build_factor_matrix(int z, double* fac) { int t; int s; *************** *** 350,359 **** } } ! static void add_crosses_to_factor_matrix(int z, float* fac, int packed_length) { ! float val; int t; int s; int i; --- 350,359 ---- } } ! static void add_crosses_to_factor_matrix(int z, double* fac, int packed_length) { ! double val; int t; int s; int i; *************** *** 368,381 **** } } ! static int pack_out_missing( float* data_in, unsigned char** missing, int tdim, int nblocks, ! int z, float* data_out ) { int iblock; int t; ! float* in; ! float* out; int n_valid; /* This program works on data for which the block index (corresponding --- 368,381 ---- } } ! static int pack_out_missing( double* data_in, unsigned char** missing, int tdim, int nblocks, ! int z, double* data_out ) { int iblock; int t; ! double* in; ! double* out; int n_valid; /* This program works on data for which the block index (corresponding *************** *** 399,411 **** return n_valid; /* note that this should be doubled for complex! */ } ! static void mean_correct_factors( float* factors, float* factor_means, int nfactors, int length ) { int ifactor; double sum; int i; ! float* runner; for (ifactor=0; ifactorchunk[mse_out->chunk_offset]= sse/(float)sse_df; mse_out->chunk_offset++; } } --- 538,549 ---- for (src_id=0; src_idchunk[mse_out->chunk_offset]= sse/(double)sse_df; mse_out->chunk_offset++; } } *************** *** 575,582 **** unsigned char** missing= NULL; int tseries_length_packed= 0; int factors_length_packed= 0; ! float mean; ! float mean_variance; int mse_flag= 0; progname= argv[0]; --- 575,582 ---- unsigned char** missing= NULL; int tseries_length_packed= 0; int factors_length_packed= 0; ! double mean; ! double mean_variance; int mse_flag= 0; progname= argv[0]; *************** *** 685,696 **** if (v_flag) { /* Do one voxel and output the results */ ! float* data; data= mri_get_chunk(Input, "images", nimages, ((((v_k*ydim) + v_j)*xdim) + v_i)*nimages, ! MRI_FLOAT); if (!data) Abort("%s: failed to read voxel %d %d %d!\n", progname,v_i,v_j,v_k); --- 685,696 ---- if (v_flag) { /* Do one voxel and output the results */ ! double* data; data= mri_get_chunk(Input, "images", nimages, ((((v_k*ydim) + v_j)*xdim) + v_i)*nimages, ! MRI_DOUBLE); if (!data) Abort("%s: failed to read voxel %d %d %d!\n", progname,v_i,v_j,v_k); *************** *** 720,726 **** } } else { ! float* data; int read_offset= 0; for (k=0; kdimstr= NULL; + f->datatype= NULL; f->length= 0; } else { safe_copy(keybuf, f->chunk); + safe_concat(keybuf, ".datatype"); + if ( !mri_has( f->file->ds, keybuf ) ) + Abort( "%s needs %s info for input file <%s>!\n", progname, keybuf, + f->file->fname ); + f->datatype= strdup(mri_get_string( f->file->ds, keybuf )); + + safe_copy(keybuf, f->chunk); safe_concat(keybuf, ".dimensions"); if ( !mri_has( f->file->ds, keybuf ) ) Abort( "%s needs %s info for input file <%s>!\n", progname, keybuf, *************** *** 360,373 **** char* default_fname_in, char* default_chunk_in, int rw) { - char keybuf[KEYBUF_SIZE]; char fname[512], chunkname[512]; - char tbuf[64]; - char* dimstr; - char* crunner; - long l; MRIChunk* f; - int created= 0; parse_fname(fname_raw, default_fname_in, fname, default_chunk_in, chunkname); f= mriChunk_open_file( mriFile_open(fname,rw), chunkname, rw ); --- 369,376 ---- *************** *** 382,388 **** char tbuf[64]; char* dimstr; char* crunner; ! long l; if (strlen(f->chunk)>KEYBUF_SIZE-64) Abort( "%s: chunk name too long!\n", progname ); --- 385,391 ---- char tbuf[64]; char* dimstr; char* crunner; ! long long l; if (strlen(f->chunk)>KEYBUF_SIZE-64) Abort( "%s: chunk name too long!\n", progname ); *************** *** 425,442 **** mriFile_unref( f->file ); } ! static void mriChunk_read(float* dest, MRIChunk* f, int n) { int nread= 0; int thisblock; ! int offset; offset= f->offset; while (n>nread) { thisblock= (n-nread > f->length - offset) ? (f->length - offset) : n-nread; if (!(mri_read_chunk(f->file->ds, f->chunk, thisblock, offset, ! MRI_FLOAT, dest+nread))) Abort("%s: bad mri_read_chunk on infile <%s>; length %d, offset %d", f->file->fname,thisblock,offset); offset += thisblock; --- 428,445 ---- mriFile_unref( f->file ); } ! static void mriChunk_read(double* dest, MRIChunk* f, int n) { int nread= 0; int thisblock; ! long long offset; offset= f->offset; while (n>nread) { thisblock= (n-nread > f->length - offset) ? (f->length - offset) : n-nread; if (!(mri_read_chunk(f->file->ds, f->chunk, thisblock, offset, ! MRI_DOUBLE, dest+nread))) Abort("%s: bad mri_read_chunk on infile <%s>; length %d, offset %d", f->file->fname,thisblock,offset); offset += thisblock; *************** *** 447,464 **** n, f->file->fname, f->chunk, (long)f->offset); } ! static void mriChunk_write(float* source, MRIChunk* f, int n) { int nwritten= 0; int thisblock; ! int offset; offset= f->offset; while (n>nwritten) { thisblock= (n-nwritten > f->length - offset) ? (f->length - offset) : n-nwritten; mri_write_chunk(f->file->ds, f->chunk, thisblock, offset, ! MRI_FLOAT, source+nwritten); offset += thisblock; if (offset>=f->length) offset= 0; nwritten += thisblock; --- 450,467 ---- n, f->file->fname, f->chunk, (long)f->offset); } ! static void mriChunk_write(double* source, MRIChunk* f, int n) { int nwritten= 0; int thisblock; ! long long offset; offset= f->offset; while (n>nwritten) { thisblock= (n-nwritten > f->length - offset) ? (f->length - offset) : n-nwritten; mri_write_chunk(f->file->ds, f->chunk, thisblock, offset, ! MRI_DOUBLE, source+nwritten); offset += thisblock; if (offset>=f->length) offset= 0; nwritten += thisblock; *************** *** 467,473 **** n, f->file->fname, f->chunk, (long)f->offset); } ! static void mriChunk_advance(MRIChunk* f, int n) { f->offset= (f->offset + n) % f->length; if (gbl_debug_flag) fprintf(stderr,"Advance %s:%s; offset now %ld\n", --- 470,476 ---- n, f->file->fname, f->chunk, (long)f->offset); } ! static void mriChunk_advance(MRIChunk* f, long long n) { f->offset= (f->offset + n) % f->length; if (gbl_debug_flag) fprintf(stderr,"Advance %s:%s; offset now %ld\n", *************** *** 568,574 **** safe_copy(keybuf,Params->chunk); safe_concat(keybuf,".datatype"); ! mri_set_string( Params->file->ds, keybuf, "float32" ); dimstr= Input->dimstr; dimstr= strchr(dimstr,'t') + 1; /* skip to past t */ --- 571,580 ---- safe_copy(keybuf,Params->chunk); safe_concat(keybuf,".datatype"); ! if (!strcmp(Input->datatype,"float64")) ! mri_set_string( Params->file->ds, keybuf, "float64" ); ! else ! mri_set_string( Params->file->ds, keybuf, "float32" ); dimstr= Input->dimstr; dimstr= strchr(dimstr,'t') + 1; /* skip to past t */ *************** *** 706,712 **** return factor; } ! static void load_factors(float* factor_raw, float* factor_scratch, int nfactors, int fac_block, MRIChunk** Factor) { int ifac; --- 712,718 ---- return factor; } ! static void load_factors(double* factor_raw, double* factor_scratch, int nfactors, int fac_block, MRIChunk** Factor) { int ifac; *************** *** 731,737 **** } ! static void apply_scale( float* data, float* scale, int complex_flag, int n, int ncols ) { int icol; --- 737,743 ---- } ! static void apply_scale( double* data, double* scale, int complex_flag, int n, int ncols ) { int icol; *************** *** 757,763 **** } } ! static void apply_unscale( float* data, float* scale, int complex_flag, int n, int ncols ) { int icol; --- 763,769 ---- } } ! static void apply_unscale( double* data, double* scale, int complex_flag, int n, int ncols ) { int icol; *************** *** 783,797 **** } } ! static int pack_out_missing( float* data_in, unsigned char** missing, int tdim, int complex_flag, int nblocks, ! int z, float* data_out ) { int iblock; int t; ! float* in; ! float* out; int n_valid; /* This program works on data for which the block index (corresponding --- 789,803 ---- } } ! static int pack_out_missing( double* data_in, unsigned char** missing, int tdim, int complex_flag, int nblocks, ! int z, double* data_out ) { int iblock; int t; ! double* in; ! double* out; int n_valid; /* This program works on data for which the block index (corresponding *************** *** 830,844 **** return n_valid; /* note that this should be doubled for complex! */ } ! static void pack_in_missing( float* data_in, unsigned char** missing, int tdim, int complex_flag, int nblocks, ! int z, float fill_val, float* data_out ) { int iblock; int t; ! float* in; ! float* out; /* This program works on data for which the block index (corresponding * to the factor) varies faster than the tdim index (corresponding to --- 836,850 ---- return n_valid; /* note that this should be doubled for complex! */ } ! static void pack_in_missing( double* data_in, unsigned char** missing, int tdim, int complex_flag, int nblocks, ! int z, double fill_val, double* data_out ) { int iblock; int t; ! double* in; ! double* out; /* This program works on data for which the block index (corresponding * to the factor) varies faster than the tdim index (corresponding to *************** *** 907,922 **** return result; } ! static void format_output_params( float* params_out, float* params, ! float* istdv_proj, ! float mean, float mean_variance, int nsamples, int nparams, int nparams_out, int nfactors ) { int i,j; int o_offset= 0; int i_offset= 0; ! float factor_scale; /* reality check, against n_output_params() */ if (nparams_out != n_output_params(nfactors)) { --- 913,928 ---- return result; } ! static void format_output_params( double* params_out, double* params, ! double* istdv_proj, ! double mean, double mean_variance, int nsamples, int nparams, int nparams_out, int nfactors ) { int i,j; int o_offset= 0; int i_offset= 0; ! double factor_scale; /* reality check, against n_output_params() */ if (nparams_out != n_output_params(nfactors)) { *************** *** 928,934 **** * freedom; there are really (nfactors+1) constraints. Fortunately * we can easily rescale. */ ! factor_scale= ((float)(nsamples-nfactors))/((float)(nsamples-(nfactors+1))); params_out[o_offset++]= mean; for (i=0; i #include ! static char rcsid[] = "$Id: libarray.c,v 1.3 1999/07/07 20:14:47 welling Exp $"; /* Alloc2DFloatArray dynamically allocates a 2-D array of floats with dimensions [d1][d2]. The actual --- 25,31 ---- #include #include ! static char rcsid[] = "$Id: libarray.c,v 1.4 2002/04/04 18:54:42 welling Exp $"; /* Alloc2DFloatArray dynamically allocates a 2-D array of floats with dimensions [d1][d2]. The actual *************** *** 129,131 **** --- 129,204 ---- free(a); } } + + /* Alloc2DDoubleArray dynamically allocates a 2-D array + of doubles with dimensions [d1][d2]. The actual + contents of the array are stored contiguously so + that block operations may be performed on the + memory chunk holding the array. */ + double ** + Alloc2DDoubleArray (int d1, + int d2) + { + int i; + double **a; + double *contents; + + contents = (double *) malloc(d1*d2*sizeof(double)); + a = (double **) malloc(d1 * sizeof(double *)); + for (i = 0; i < d1; ++i) + a[i] = contents + i*d2; + return(a); + } + + /* Free2DDoubleArray deallocates a 2-D array created + by Alloc2DDoubleArray */ + void Free2DDoubleArray (double **a) + { + if (a != NULL) + { + free(a[0]); + free(a); + } + } + + /* Alloc3DDoubleArray dynamically allocates a 3-D array + of doubles with dimensions [d1][d2][d3]. The actual + contents of the array are stored contiguously so + that block operations may be performed on the + memory chunk holding the array. */ + double *** + Alloc3DDoubleArray (int d1, + int d2, + int d3) + { + int i, j; + double ***a; + double **index; + double *contents, *offset; + + contents = (double *) malloc(d1*d2*d3*sizeof(double)); + index = (double **) malloc(d1*d2*sizeof(double *)); + a = (double ***) malloc(d1*sizeof(double **)); + for (i = 0; i < d1; ++i) + { + a[i] = index + i*d2; + offset = contents + i*d2*d3; + for (j = 0; j < d2; ++j) + a[i][j] = offset + j*d3; + } + return(a); + } + + /* Free3DDoubleArray deallocates a 3-D array created + by Alloc3DDoubleArray */ + void Free3DDoubleArray (double ***a) + { + if (a != NULL) + { + free(a[0][0]); + free(a[0]); + free(a); + } + } + +