#!/usr/bin/env AFNI_Batch_R #iplots are quite cool, # but when it comes to plotting multiple time series (say 576 x 42) # it gets reaaaaal slow. JAVA SUCKS. #tkrplot has the stupid feature of not allowing autoresize of the plot. # you can do it manually by setting hscale and vscale, but I'd rather # play with a copperhead instead #Clean up rm(list = ls()) first.in.path <- function(file) { ff <- paste(strsplit(Sys.getenv('PATH'),':')[[1]],'/', file, sep='') ff<-ff[lapply(ff,file.exists)==TRUE]; #cat('Using ', ff[1],'\n'); return(gsub('//','/',ff[1], fixed=TRUE)) } source(first.in.path('AFNIio.R')) source(first.in.path('AFNIplot.R')) ExecName <- 'ExamineXmat' greeting.ExamineXmat <- function () return( "#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ================== Welcome to ExamineXmat.R ================== A program to examine the design matrix #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" ) reference.ExamineXmat <- function () return( " Ziad S. Saad (SSCC/NIMH/NIH)" ) #The help function for ExamineXmat batch (command line mode) help.ExamineXmat.opts <- function (params, alpha = TRUE, itspace=' ', adieu=FALSE) { intro <- ' Usage: ------ ExamineXmat is a program for examining the design matrix generated by 3dDeconvolve. The examination can be done interactively, by entering your selections in a GUI (-interactive). Alternately, you can send the output to an image file (-prefix). The title of the plot in the image contains the filename and the number of regressors shown/total number. The subtitle contains the various condition numbers of interest: Rall for the entire matrix Rall-motion for the matrix without motion parameters Rall-roni for the matrix without any regressors of no interest Rviewed for the part of the matrix displayed in the graph More help is available in interactive usage mode. ' ex1 <- " Examples --- : -------------- #Get some sample data curl -o demo.X.xmat.1D \\ afni.nimh.nih.gov/pub/dist/edu/data/samples/X.xmat.1D #PDF output may not work on some machines ExamineXmat -prefix t1.pdf -input demo.X.xmat.1D ExamineXmat -prefix t2.jpg -input demo.X.xmat.1D -select ALL_TASKS ExamineXmat -prefix t3.png -input demo.X.xmat.1D -select tneg tpos ExamineXmat -prefix t4.jpg -input demo.X.xmat.1D -select MOTION tneg ExamineXmat -prefix t5.jpg -input demo.X.xmat.1D -select tneg 3 35:38 #interactive mode ExamineXmat ExamineXmat -input demo.X.xmat.1D ExamineXmat -input demo.X.xmat.1D -select tneg #To save the last image you see interactively ExamineXmat -input demo.X.xmat.1D -interactive \\ -prefix t6.jpg -select tneg tpos " parnames <- names(params) ss <- vector('character') if (alpha) { parnames <- sort(parnames) ss <- paste('Options in alphabetical order:\n', '------------------------------\n', sep='') } else { ss <- paste('Options:\n', '--------\n', sep='') } for (ii in 1:length(parnames)) { op <- params[parnames[ii]][[1]] if (!is.null(op$help)) { ss <- c(ss , paste(itspace, op$help, sep='')); } else { ss <- c(ss, paste(itspace, parnames[ii], '(no help available)\n', sep='')); } } ss <- paste(ss, sep='\n'); cat(intro, ex1, ss, reference.ExamineXmat(), sep='\n'); if (adieu) exit.AFNI(); } init.ExamineXmat.lop <- function () { lop <- list() lop$input = NULL; lop$interactive = FALSE; lop$pprefix=NULL; lop$cprefix=NULL; lop$select = 'ALL_TASKS'; lop$verb=0 return(lop) } set.both.prefix <- function (lop, pp) { pp <- parse.name(pp) lop$pprefix <- paste(pp$path,pp$name_noext, ".plt",pp$ext, sep='') lop$cprefix <- paste(pp$path,pp$name_noext, ".cor",pp$ext, sep='') return(lop) } #Change command line arguments into an options list read.ExamineXmat.opts.batch <- function (args=NULL, verb = 0) { params <- list ( '-input' = apl(n = c(1), d = NA, h = paste( "-input 1Dfile: xmat file to plot\n" ) ), '-select' = apl(n = c(1, Inf), d = 'ALL_TASKS', h = paste( "-select SELECTION_STRING: What to plot.\n", " Selection Strings:\n", pad.string.lines(ExamineXmat.help.select(),' '), "\n Alternately you can specify special strings:\n", " 'ALL': The entire matrix\n", " 'ALL_TASKS': All task regressors\n", " 'RONI': All regressors of no interest (baseline+motion)\n", " 'BASE': All baseline regressors\n", " 'MOT': All motion regressors\n" ) ), '-prefix' = apl(n = 1, d = NA, h = paste( "-prefix PREFIX: Prefix of plot image and cor image \n" ) ), '-pprefix' = apl(n = 1, d = NA, h = paste( "-pprefix PPREFIX: Prefix of plot image only \n" ) ), '-cprefix' = apl(n = 1, d = NA, h = paste( "-cprefix CPREFIX: Prefix of cor image only \n" ) ), '-interactive' = apl(n = 0, d = NA, h = paste( "-interactive: Run examine Xmat in interactive mode\n", " This is the default if -prefix is not given.\n", " If -interactive is used with -prefix, the last\n", " plot you see is the plot saved to file.\n" ) ), '-verb' = apl(n=1, d = 0, h = paste( "-verb VERB: VERB is an integer specifying verbosity level.\n", " 0 for quiet (Default). 1 or more: talkative.\n" ) ), '-help' = apl(n=0, h = '-help: this help message\n'), '-h' = apl(n=0, h = '-h: this help message\n'), '-show_allowed_options' = apl(n=0, h= "-show_allowed_options: list of allowed options\n" ), '-msg.trace' = apl(n=0, h= "-msg.trace: Output trace information along with errors and notices\n" ) ); ops <- parse.AFNI.args(args, params, other_ok=FALSE, verb=verb ) if (verb) show.AFNI.args(ops, verb=verb-1, hstr=''); if (is.null(ops)) { errex.AFNI('Error parsing arguments. See ExamineXmat -help for details.'); } #Parse dems options #initialize with defaults lop <- init.ExamineXmat.lop() #Get user's input for (i in 1:length(ops)) { opname <- strsplit(names(ops)[i],'^-')[[1]]; opname <- opname[length(opname)]; switch(opname, input = lop$input <- ops[[i]], prefix = lop <- set.both.prefix(lop,ops[[i]]), pprefix = lop$pprefix <- ops[[i]], cprefix = lop$cprefix <- ops[[i]], interactive = lop$interactive <- TRUE, select = lop$select <- ops[[i]], verb = lop$verb <- ops[[i]], help = help.ExamineXmat.opts(params, adieu=TRUE), h = help.ExamineXmat.opts(params, adieu=TRUE), show_allowed_options = show.AFNI.args(ops, verb=0, hstr="ExamineXmat's",adieu=TRUE), msg.trace = set.AFNI.msg.trace(TRUE) ) } return(lop) }# end of read.ExamineXmat.opts.batch #Change options list to ExamineXmat variable list process.ExamineXmat.opts <- function (lop, verb = 0) { if (is.null(lop$pprefix)) lop$interactive <- TRUE return(lop) } imgdev <<- NULL plotdev <<- NULL #The plotting function show_xmat <- function (xmat, isel=1:1:ncol(xmat), descr="", interactive=TRUE) { #cat ('isel in show_xmat') #show(isel) if (length(isel) == 0) { err.AFNI("Nothing to show"); return(0) } #Set colors based on group and offset by 2 so that -1 and 0 get colored colvec <- cg <- attr(xmat,'ColumnGroups') #first color regressors of no interest colvec[cg == -1] = 1 colvec[cg == 0 ] = 2 N_basegroup = 2 #2 types of baseline (RONI) #now color regressors of interest, skip yellow collist = c(3,4,5,6,8) #all but 1st two and yellow N_collist = length(collist) for (i in min(colvec[colvec > 0]):max(colvec[colvec > 0])) { colvec[cg == i] <- collist[(i-N_basegroup) %% N_collist + 1] } if (length(isel) > 10) { text.lym <- colnames(xmat) text.rym <- paste('c', seq(0,ncol(xmat)-1,1), sep='') } else { text.lym <- NULL text.rym <- paste('c', seq(0,ncol(xmat)-1,1), sep='') } descr <- paste( 'Xmat ', attr(xmat,'FileName'), '\n', 'Baseline indices: ', xmat.base.index(xmat)-1, '\n', 'Motion indices: ', xmat.motion.index(xmat)-1, '\n', 'Task indices: ', xmat.alltasks.index(xmat)-1, '\n'); if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) }else{ view_cond <- 0 } if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- 0 } stit = paste ( 'Rall: ', sprintf('%.2f', all_cond), '', 'Rall-motion: ', sprintf('%.2f', all_cond_no_mot), '', 'Rall-roni: ', sprintf('%.2f', stim_cond),'', 'Rviewed: ', sprintf('%.2f', view_cond),''); maintit <- sprintf('%s\n(%d/%d)', trim.string(attr(xmat,'FileName'), 48), length(isel), ncol(xmat)) maintitcorr <- sprintf('cor(%s)\n(%d/%d)', trim.string(attr(xmat,'FileName'), 48), length(isel), ncol(xmat)) if (!is.null(plotdev) && plotdev != 0) { dev.set(plotdev) } plotdev <<- plot.1D( dmat = xmat, dmat.err=NULL, dmat.colsel = isel, col.ystack = TRUE, ttl.main = maintit, ttl.sub=stit, multi.ncol=1, col.mean.line = FALSE, col.color = colvec, col.plot.char = -1, col.plot.type = 'l', leg.show = FALSE, col.text.lym = text.lym, col.text.lym.at = 'YOFF', col.text.rym = text.rym, col.text.rym.at = 'YOFF', xax.label = 'Time', prefix = lop$pprefix, nodisp = !interactive ) #update report window, if interactive mode if (interactive) { #title (paste('Xmat: ', attr(xmat,'FileName'),'\n')); title (''); ttc <<- condition_report(ttc, xmat, isel, descr=descr) } #Need a device for correlation matrix display imgdev <<- image.corr.1D(dmat = xmat, dmat.colsel = isel, prefix = lop$cprefix, nodisp = !interactive, ttl.main = maintitcorr, dev.this = imgdev, dev.new=TRUE) return(0) } show_xmat_old <- function (xmat, isel=1:1:ncol(xmat), descr="") { #cat ('isel in show_xmat') #show(isel) #Set colors based on group and offset by 2 so that -1 and 0 get colored colvec <- cg <- attr(xmat,'ColumnGroups') #first color regressors of no interest colvec[cg == -1] = 1 colvec[cg == 0 ] = 2 N_basegroup = 2 #2 types of baseline (RONI) #now color regressors of interest, skip yellow collist = c(3,4,5,6,8) #all but 1st two and yellow N_collist = length(collist) for (i in min(colvec[colvec > 0]):max(colvec[colvec > 0])) { colvec[cg == i] <- collist[(i-N_basegroup) %% N_collist + 1] } if (length(isel) > 10) { tp = 'single' #Calculate normalized version for full display ffr <- xmat[,isel] ra = matrix(nrow=ncol(ffr), ncol=3) for (i in 1:ncol(ffr)) { ra[i,] <- as.vector(quantile(ffr[,i], c(0,1,0.5))) } offs <-vector(length = ncol(ffr), mode="numeric") for (i in 1:ncol(ffr)) { ffr[,i] <- (ffr[,i]-ra[i,1])/(ra[i,2]-ra[i,1])+1.4*i offs[i] = mean(ffr[,i]) } if (is.null(dev.list())) x11() plot(ffr, plot.type = tp, xy.labels = colnames(xmat[,isel]), xy.lines = TRUE, panel = lines, nc = 3, yax.flip = FALSE, axes = TRUE, col = colvec[isel], main = '') xoff = ((1:length(isel))%%5)*(length(xmat[,1])/frequency(xmat)/5) text (xoff, offs, colnames(xmat[,isel]), col = colvec[isel], adj=c(0,0)) } else { tp = 'multiple' if (is.null(dev.list())) x11() if (length(isel) == 1) { colm <- colvec[isel] } else { colm <- 1 #ARRRRGH only one col is used for multiple line drawings! } plot(xmat[,isel], plot.type = tp, xy.labels = colnames(xmat[,isel]), xy.lines = TRUE, panel = lines, nc = 1, yax.flip = FALSE, axes = TRUE, col = colm, main = '') } cat ( 'Xmat ', attr(xmat,'FileName'), '\n', 'Baseline indices: ', xmat.base.index(xmat)-1, '\n', 'Motion indices: ', xmat.motion.index(xmat)-1, '\n', 'Task indices: ', xmat.alltasks.index(xmat)-1, '\n'); if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) }else{ view_cond <- 0 } if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- 0 } stit = paste ( 'Rall : ', sprintf('%.2f', all_cond), '\n', 'Rall-motion : ', sprintf('%.2f', all_cond_no_mot), '\n', 'Rall-roni : ', sprintf('%.2f', stim_cond),'\n', 'Rviewed : ', sprintf('%.2f', view_cond),'\n'); #title (paste('Xmat: ', attr(xmat,'FileName'),'\n')); title (''); #update report window ttc <<- condition_report(ttc, xmat, isel, descr=descr) return(0) } text_window_kill <- function() { tkdestroy(ttc) ttc <<-NA } condition_report <- function(tt=NA, xmat, isel=1:1:ncol(xmat), descr="") { if (is.na(tt)[1]) { #cat ('Have no widget\n') tt <- tktoplevel() #Catch the destroy button to reset ttc to NA tcl( "wm", "protocol", tt, "WM_DELETE_WINDOW", function()text_window_kill()) scr <- tkscrollbar(tt, repeatinterval=5, command=function(...)tkyview(txt,...)) #make sure txt is being kept for next call txt <<- tktext( tt,bg="white",font="courier", yscrollcommand=function(...)tkset(scr,...)) tkgrid(txt,scr) tkgrid.configure(scr,sticky="ns") #Numbers that will stay constant if (length(xmat.alltasks.index(xmat))) { stim_cond <- kappa(xmat[,xmat.alltasks.index(xmat)], exact=TRUE) } else { stim_cond <- 0 } all_cond <- kappa(xmat, exact=TRUE) if (length(xmat.motion.index(xmat))) { all_cond_no_mot <- kappa(xmat[,-xmat.motion.index(xmat)], exact=TRUE) } else { all_cond_no_mot <- all_cond } #attr(xmat,'TR') #attr(xmat,'dim') stitcom = paste ( 'Matrix file: ', attr(xmat,'FileName'), '\n', '* TR: ', attr(xmat,'TR'), '\n', '* N: ', length(xmat[,1]), '\n', '* Baseline indices: ', paste(xmat.base.index(xmat)-1, collapse=", "), '\n', '* Motion indices: ', paste(xmat.motion.index(xmat)-1, collapse=", "), '\n', '* Task indices: ', paste(xmat.alltasks.index(xmat)-1, collapse=", "), '\n', '* Task Labels: ', paste( attr(xmat,'TaskNames'), collapse = ", "),'\n', '* Standard condition numbers:', '\n', ' Rall : ', sprintf('%.2f', all_cond), '\n', ' Rall-motion : ', sprintf('%.2f', all_cond_no_mot), '\n', ' Rall-roni : ', sprintf('%.2f', stim_cond),'\n', '=============================================\n'); } else { #cat ('reusing widget\n') stitcom = paste ('User Selection:', paste(descr, collapse=" ") , '\n') } #cat ('Calculating conditions\n') if (length(isel)) { view_cond <- kappa(xmat[,isel], exact=TRUE) } else { view_cond = 0 } stit = paste ( stitcom, ' Rviewed : ', sprintf('%.2f', view_cond),'\n', '----------------------------------------------\n'); #cat ('inserting string\n') tkconfigure(txt, state="normal") #allow editing tkinsert(txt,"end",stit) tkconfigure(txt, state="disabled") #no more editing tkfocus(txt) return(tt) } sel_gui <- function (tt, xmat) { require(tcltk) Name <- tclVar("") entry.Name <-tkentry(tt,width="20",textvariable=Name) tkgrid(tklabel(tt, text=paste( "Enter indices or label identifiers\n", "example: 0:3, ", attr(xmat,'TaskNames')[1],", 7 9") ) ) tkgrid(entry.Name) OnOK <- function() { #cat ('OnOK\n') NameVal <- gsub("[,;\'\"]", " ",tclvalue(Name)) show_sel(NameVal,xmat) tkfocus(tt) } OnDone <- function() { tkdestroy(tt) exit_flag <<- TRUE } OnHelp <- function() { ExamineXmat.help() } OK.but <-tkbutton(tt,text=" OK ",command=OnOK) Done.but <- tkbutton(tt,text=" Done ",command=OnDone) Help.but <- tkbutton(tt,text=" Help ",command=OnHelp) tkbind(entry.Name, "", OnOK) tkgrid(OK.but, Done.but, Help.but) tkfocus(tt) exit_flag <<- FALSE } ExamineXmat.help.select <- function() { sh <- paste ( 'To select regressors, you can use regressor indices or regressor \n', 'labels. For example, say your tasks (and therefore regressors) \n', 'are labeled "house", "face", "random"\n', 'Then to view house and face regressors only, select:\n', '"house, face" or "house face" or "h fa" etc.\n', 'You can also use regressor indices (start at 0)\n', 'So you can select "0, 5, 10:12" for regressors 0, 5, 10, 11, and 12\n', 'You can also combine strings and integers in the same selection.\n', 'Commas, semicolons, and quotes are ignored.\n' ) return(sh) } ExamineXmat.help <- function() { sh <- paste ( '===============================================================\n', 'ExamineXmat Help:\n', '-----------------\n', 'ExamineXmat.R is an interactive tool to examine a design matrix.\n', 'The interface consists of 3 windows: Graph, Control, Results.\n', '\n', 'The Graph window displays the entire design matrix or a subset \n', 'designated by the user.\n', 'The Control window allows users to designate which regressors\n', 'to display.\n', 'The Results window displays properties of the design matrix in \n', 'its entirety and the condition number of the matrix formed \n', 'by the designated (or viewed in the graph) regressors only.\n', '\n', ExamineXmat.help.select(), '\n', 'When viewing more than 10 regressors, the regressors get automatically\n', 'rescaled and displayed on the same axis. In this case, the amplitude\n', 'of the Y axis is of no relevance. In this display mode, baseline \n', 'regressors are black, motion are red, and each task group gets a distinct\n', 'color.\n', 'When 10 or less regressors are displayed simultaneously, each one\n', 'gets its own subplot and the y axis represents the regressor amplitude\n', 'as it is in the design matrix. Also, in this display mode, all regressors\n', 'are drawn in black.\n', '===============================================================\n', '\n') cat (sh) #tkmessageBox(message=sh) #ugly... #use reults output instead tkconfigure(txt, state="normal") #allow editing tkinsert(txt,"0.0",sh) tkconfigure(txt, state="disabled") #no more editing tkfocus(txt) #tkmessageBox(message="See text in Results window") } show_sel <- function(sel, xmat) { ilst <- xmat.select.indices(sel, xmat) if (length(ilst)) { ilst = unique(ilst[ilst >0 & ilst <=ncol(xmat)]) if (length(ilst)) { cat ('selection is: ', ilst-1, '\n') show_xmat(xmat, ilst, descr=sel) } } } main_loop <- function (fn, select=NULL) { #read file if (is.null(lop$xmat <- read.AFNI.xmat(fn))) { err.AFNI("Failed to load xmat"); return(1); } #top level widget ttc <<- NA #Show all of lop$xmatr isel <- NULL if (!is.null(select)) { isel <- xmat.select.indices(select, lop$xmat) if (length(isel)==0) { warn.AFNI(sprintf("Nothing found for selection %s. Showing all.", select)) isel<-NULL; } } if (is.null(isel)) isel <- 1:1:ncol(lop$xmat) show_xmat(lop$xmat, isel) #begin the selection toy, this GUI controls the quitting etc. sel_gui(tktoplevel(), lop$xmat) #dunno what to do for event loop, all I can do is sleep #Without this, the batch mode exits immediately while (!exit_flag) { Sys.sleep(0.5) } return(0) } ################################################################################# ####################### Begin ExamineXmat main ################################## ################################################################################# if (!exists('.DBG_args')) { args = (commandArgs(TRUE)) rfile <- first.in.path(sprintf('%s.R',ExecName)) save(args, rfile, file=".ExamineXmat.dbg.AFNI.args", ascii = TRUE) } else { note.AFNI("Using .DBG_args resident in workspace"); args <- .DBG_args } if (!length(args)) { BATCH_MODE <<- 0 err.AFNI("No parameters"); } else { if (!exists('.DBG_args')) { BATCH_MODE <<- 1 } else { BATCH_MODE <<- 0 } if (is.null(lop <- read.ExamineXmat.opts.batch(args, verb = 0))) { stop('Error parsing input'); } if (is.null(lop <- process.ExamineXmat.opts(lop, verb = lop$verb))) { stop('Error processing input'); } } if (lop$verb) { str(lop); } if (lop$interactive) { require('graphics') require('tcltk') if (is.null(lop$input)) { lop$input <- tclvalue( tkgetOpenFile( filetypes = "{{Xmat Files} {.xmat.1D}} {{All files} *}", title = 'Choose design matrix file')) } if (nchar(lop$input)) { if (main_loop(lop$input, lop$select)) cat ('error') } } else { if (is.null(lop$xmat <- read.AFNI.xmat(lop$input))) { errex.AFNI("Failed to load xmat"); } isel <- NULL if (!is.null(lop$select)) { isel <- xmat.select.indices(lop$select, lop$xmat) if (length(isel)==0) { errex.AFNI(sprintf("Nothing found for selection %s.", lop$select)) isel<-NULL; } } if (is.null(isel)) { isel <- 1:1:ncol(lop$xmat) } show_xmat(lop$xmat, isel, interactive=FALSE) }