1#!/usr/bin/env AFNI_Batch_R 2 3#Clean up 4rm(list = ls()) 5 6if (1) { 7 #MASS's fitdistr, does not converge always 8 #library(MASS) 9 #fitdistrplus may be better 10 #install.packages("fitdistrplus") 11 if (!require(fitdistrplus)) { 12 require(MASS) 13 } 14} 15 16first.in.path <- function(file) { 17 ff <- paste(strsplit(Sys.getenv('PATH'),':')[[1]],'/', file, sep='') 18 ff<-ff[lapply(ff,file.exists)==TRUE]; 19 #cat('Using ', ff[1],'\n'); 20 return(gsub('//','/',ff[1], fixed=TRUE)) 21} 22source(first.in.path('AFNIio.R')) 23source(first.in.path('Signatures.R')) 24source(first.in.path('AFNIplot.R')) 25 26#And load R's libsvm 27libLoad("e1071") 28 29 30ExecName <- '3dSignatures' 31 32greeting.SigsClassify <- function () 33 return( "#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 34 ================== Welcome to SigsClassify.R ================== #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" 35 ) 36 37reference.SigsClassify <- function () 38 return( 39"" 40 ) 41 42 43#The help function for SigsClassify batch (command line mode) 44help.SigsClassify.opts <- function (params, alpha = TRUE, itspace=' ', adieu=FALSE) { 45 46 intro <- 47' 48Usage: 49------ 50 3dSignatures is a program creating voxel membership priors based on 51 a collection of features termed "signatures".' 52 53 ex1 <- 54" 55Example 1 --- : 56-------------------------------- 57 SigsClassify -input IBSR_01.SubClasses.skew.1D.repsig 58 59 SigsClassify -input IBSR_01.SubClasses.skew.1D.repsig \\ 60 -ColumnGroups 'R:c(rep(1,10), rep(2,10), rep(3,10))' \\ 61 -GroupLabels CSF GM WM \\ 62 -NoZeros -Nocol.ystack \\ 63 -trainsuffix ZePlotIzSaved.pdf 64" 65 66 parnames <- names(params) 67 ss <- vector('character') 68 if (alpha) { 69 parnames <- sort(parnames) 70 ss <- paste('Options in alphabetical order:\n', 71 '------------------------------\n', sep='') 72 } else { 73 ss <- paste('Options:\n', 74 '--------\n', sep='') 75 } 76 for (ii in 1:length(parnames)) { 77 op <- params[parnames[ii]][[1]] 78 if (!is.null(op$help)) { 79 ss <- c(ss , paste(itspace, op$help, sep='')); 80 } else { 81 ss <- c(ss, paste(itspace, parnames[ii], 82 '(no help available)\n', sep='')); 83 } 84 } 85 ss <- paste(ss, sep='\n'); 86 cat(intro, ex1, ss, reference.SigsClassify(), sep='\n'); 87 88 if (adieu) exit.AFNI(); 89} 90 91parse.SigsClassify.ColumnGroups <- function(op) { 92 mm <- read.AFNI.matrix(op) 93 return(mm) 94} 95 96init.SigsClassify.lop <- function () { 97 lop <- list() 98 lop$tset <- list() 99 lop$Tset <- list() 100 lop$odir = NULL; 101 lop$opref = 'test'; 102 lop$volsuffix = NULL; 103 lop$trainsuffix = 'train'; 104 lop$testsuffix = 'test'; 105 lop$thiskeyset= NULL; 106 lop$thislabset= NULL; 107 lop$labeltablefile = NULL; 108 lop$groupedlabeltablefile = NULL; 109 lop$verb=0; 110 lop$no_X11=FALSE; 111 lop$doprob = TRUE; 112 lop$no_tune=FALSE; 113 lop$svmscale=TRUE; 114 lop$reuse_train=FALSE; 115 lop$load_train=NULL; 116 lop$samples_frac <- NULL; 117 lop$max_samples <- NULL; 118 lop$min_samples <- NULL; 119 lop$method <- 'svm' 120 lop$uidp <- '' 121 return(lop) 122} 123 124parse.TeSetSF.opts <- function (iset, op) { 125 iset <- c(iset, 126 list( 127 list(subj=op[1], 128 anat=NULL, 129 anatname=NULL, 130 sig=NULL, 131 signame=op[2], 132 lab=NULL, 133 labname=NULL, 134 labselecs=NULL))) 135 return(iset) 136} 137parse.TeSetSFL.opts <- function (iset, op) { 138 iset <- c(iset, 139 list( 140 list(subj=op[1], 141 anat=NULL, 142 anatname=NULL, 143 sig=NULL, 144 signame=op[2], 145 lab=NULL, 146 labname=op[3], 147 labselecs=NULL))) 148 return(iset) 149} 150 151parse.TeSetSAFL.opts <- function (iset, op) { 152 iset <- c(iset, 153 list( 154 list(subj=op[1], 155 anat=NULL, 156 anatname=op[2], 157 sig=NULL, 158 signame=op[3], 159 lab=NULL, 160 labname=op[4], 161 labselecs=NULL))) 162 return(iset) 163} 164 165parse.TrSetSFL.opts <- function (iset, op) { 166 iset <- c(iset, 167 list( 168 list(subj=op[1], 169 sig=NULL, 170 signame=op[2], 171 lab=NULL, 172 labname=op[3], 173 labselecs=NULL))) 174 return(iset) 175} 176 177parse.CSet.opts <- function (iset, op) { 178 iset <- c(iset, 179 list( 180 list(subj=op[1], 181 sig=NULL, 182 signame=op[2], 183 lab=NULL, 184 labname=NULL, 185 labselecs=NULL))) 186 return(iset) 187} 188 189#Change command line arguments into an options list 190read.SigsClassify.opts.batch <- function (args=NULL, verb = 0) { 191 192 params <- list ( 193 '-method' = apl(n=1, d="svm", dup=FALSE, h=paste( 194 "-method METH: Choose approach. Either 'svm', or 'prob'" 195 ) ), 196 '-train' = apl(n = 3, d = NA, dup = TRUE, h = paste( 197 "See -train_SFL" 198 ) ), 199 '-train_SFL' = apl(n = 3, d = NA, dup = TRUE, h = paste( 200 "-train_SFL SUBJ FEATURES LABELS: Specify a training set's \n", 201 " subject ID, volume of signatures, and\n", 202 " volume of labels" 203 ) ), 204 '-train_suffix' = apl(n = 1, d = NA, h = paste( 205 "-train_suffix TRAINSUFFIX: Suffix of file containing training results\n" 206 ) ), 207 '-no_tune' = apl(n=0, d = 0, h = paste( 208 "-no_tune : No tuning of training model\n" 209 ) ), 210 '-test' = apl(n = 2, d = NA, dup = TRUE, h = paste( 211 "See -test_SF\n" 212 ) ), 213 '-test_SF' = apl(n = 2, d = NA, dup = TRUE, h = paste( 214 "-test_SF SUBJ FEATURES : Specify a test set's \n", 215 " subject ID, volume of signatures\n" 216 ) ), 217 '-test_SFL' = apl(n = 3, d = NA, dup = TRUE, h = paste( 218 "-test_SFL SUBJ FEATURES LABELS: Specify a test set's \n", 219 " subject ID, volume of signatures, and\n", 220 " volume of labels" 221 ) ), 222 '-test_SAFL' = apl(n = 4, d = NA, dup = TRUE, h = paste( 223 "-test_SAFL SUBJ ANAT FEATURES LABELS: Specify a test set's \n", 224 " subject ID, Anatomical, \n", 225 " volume of signatures, and volume of labels." 226 ) ), 227 '-class' = apl(n = 2, d = NA, dup = TRUE, h = paste( 228 "-test SUBJ SIGS LABELS: Specify a classification set's \n", 229 " subject ID, and volume of signatures" 230 ) ), 231 '-odir' = apl(n = 1, d = NA, dup = FALSE, h = paste( 232 "-odir ODIR: Specify a directory for output of classification results" 233 ) ), 234 '-class_labels' = apl(n = c(1,Inf), d = NA, h = paste( 235 "-class_labels CS1 CS2 ...: \n" 236 ) ), 237 '-class_keys' = apl(n = c(1,Inf), d = NA, h = paste( 238 "-class_keys Ci1 Ci2 ...: \n" 239 ) ), 240 241 '-parent_grouping' = apl(n = c(1,Inf), d = NA, h = paste( 242 "-parent_grouping ...: \n" 243 ) ), 244 245 '-parent_group_labeling' = apl(n = c(1,Inf), d = NA, h = paste( 246 "-parent_group_labeling ...: \n" 247 ) ), 248 249 '-labeltable' = apl(n = 1, d = NA, h = paste( 250 "-labeltable labletable: \n" 251 ) ), 252 '-grouped_labeltable' = apl(n = 1, d = NA, h = paste( 253 "-grouped_labeltable grouped_labletable: \n" 254 ) ), 255 '-vol_suffix' = apl(n = 1, d = NA, h = paste( 256 "-vol_suffix PREFIX: Output suffix\n" 257 ) ), 258 '-load_train' = apl(n = 1, d = NA, h = paste( 259 "-load_train PREFIX: Output suffix\n" 260 ) ), 261 '-no_X11' = apl(n=0, d = 0, h = paste( 262 "-no_X11 : No X11 display possible\n" 263 ) ), 264 '-no_prob' = apl(n=0, d = 0, h = paste( 265 "-no_prob : No probability output\n" 266 ) ), 267 '-reuse_train' = apl(n=0, d = 0, h = paste( 268 "-reuse_train : Reuse training result, if it exists\n" 269 ) ), 270 '-featscale' = apl(n=1, d = 1, h = paste( 271 "-featscale : 0/[1] scaling of features\n" 272 ) ), 273 '-min_samples' = apl(n=1, d = 100, h=paste( 274 "-min_samples MINL : Min number of training samples\n" 275 ) ), 276 '-max_samples' = apl(n=1, d = 100, h=paste( 277 "-max_samples MAX : Max number of training samples\n" 278 ) ), 279 '-sampling_fraction' = apl(n=1, d = NA, h = paste( 280 "-sampling_fraction: Fraction of voxel from each training volume\n", 281 " to use for model building.\n" 282 ) ), 283 '-uidp' = apl(n=1, d = NA, h = paste( 284 "-uidp: Unique ID prefix\n" 285 ) ), 286 '-verb' = apl(n=1, d = 0, h = paste( 287 "-verb VERB: VERB is an integer specifying verbosity level.\n", 288 " 0 for quiet (Default). 1 or more: talkative.\n" 289 ) ), 290 '-help' = apl(n=0, h = '-help: this help message\n'), 291 '-show_allowed_options' = apl(n=0, h= 292 "-show_allowed_options: list of allowed options\n" ), 293 '-msg.trace' = apl(n=0, h= 294 "-msg.trace: Output trace information along with errors and notices\n" ) 295 296 ); 297 298 ops <- parse.AFNI.args(args, params, 299 other_ok=FALSE, verb=verb 300 ) 301 if (verb) show.AFNI.args(ops, verb=verb-1, hstr=''); 302 if (is.null(ops)) { 303 errex.AFNI(paste( 304 'Error parsing arguments. See ',ExecName, 305 'SigsClassify -help for details.', sep='')); 306 } 307 308 #Parse dems options 309 #initialize with defaults 310 lop <- init.SigsClassify.lop() 311 312 #Get user's input 313 for (i in 1:length(ops)) { 314 opname <- strsplit(names(ops)[i],'^-')[[1]]; 315 opname <- opname[length(opname)]; 316 switch(opname, 317 method = lop$method <- ops[[i]], 318 train = lop$tset <- parse.TrSetSFL.opts(lop$tset, ops[[i]]), 319 train_SFL = lop$tset <- parse.TrSetSFL.opts(lop$tset, ops[[i]]), 320 test = lop$Tset <- parse.TeSetSF.opts(lop$Tset, ops[[i]]), 321 test_SF = lop$Tset <- parse.TeSetSF.opts(lop$Tset, ops[[i]]), 322 test_SFL = lop$Tset <- parse.TeSetSFL.opts(lop$Tset, ops[[i]]), 323 test_SAFL = lop$Tset <- parse.TeSetSAFL.opts(lop$Tset, ops[[i]]), 324 class = lop$Tset <- parse.CSet.opts(lop$Tset, ops[[i]]), 325 odir = lop$odir <- ops[[i]], 326 opref = lop$opref <- ops[[i]], 327 vol_suffix = lop$volsuffix <- ops[[i]], 328 train_suffix = lop$trainsuffix <- 329 strip.extension(ops[[i]], 330 c(".Rdat", ".niml.td"))$name_noext, 331 test_suffix = lop$testsuffix <- ops[[i]], 332 load_train = lop$load_train <- ops[[i]], 333 class_labels = lop$thislabset <- ops[[i]], 334 class_keys = lop$thiskeyset <- read.AFNI.matrix(ops[[i]]), 335 parent_grouping = lop$ColumnGroups <- read.AFNI.matrix(ops[[i]]), 336 parent_group_labeling = lop$GroupLabels <- ops[[i]], 337 labeltable = lop$labeltablefile <- ops[[i]], 338 grouped_labeltable = lop$groupedlabeltablefile <- ops[[i]], 339 verb = lop$verb <- ops[[i]], 340 help = help.SigsClassify.opts(params, adieu=TRUE), 341 show_allowed_options = show.AFNI.args(ops, verb=0, 342 hstr="SigsClassify's",adieu=TRUE), 343 msg.trace = set.AFNI.msg.trace(TRUE), 344 no_prob = lop$doprob <- FALSE, 345 no_X11 = lop$no_X11 <- TRUE, 346 no_tune = lop$no_tune <- TRUE, 347 reuse_train = lop$reuse_train <- TRUE, 348 sampling_fraction = lop$samples_frac <- ops[[i]], 349 min_samples = lop$min_samples <- as.numeric(ops[[i]]), 350 max_samples = lop$max_samples <- as.numeric(ops[[i]]), 351 uidp = lop$uidp <- ops[[i]], 352 featscale = lop$svmscale <- as.logical(ops[[i]]) 353 ) 354 } 355 if (lop$svmscale != FALSE && lop$svmscale != TRUE) { 356 errex.AFNI("Bad -featscale value"); 357 } 358 if (lop$no_X11 == FALSE) { 359 #Make sure that is possible or turn option off. 360 if (!length(kk <- capabilities(what='X11')) || kk[1] == FALSE) { 361 lop$no_X11 <- TRUE 362 } 363 } 364 365 if (!is.null(lop$odir) && !file_test('-d', lop$odir)) { 366 errex.AFNI(sprintf('Output directory "%s" does not exist', lop$odir)); 367 } 368 369 if (!is.null(lop$labeltablefile) && !file_test('-f', lop$labeltablefile)) { 370 errex.AFNI(sprintf('labeltable file "%s" does not exist', 371 lop$labeltablefile)); 372 } 373 374 if (!is.null(lop$thislabset)) { 375 if (!is.null(lop$thiskeyset)) { 376 errex.AFNI("You have used both -class_keys and -class_labels") 377 } 378 lop$thiskeyset <- labelkey.labeltable(lop$thislabset, lop$labeltablefile) 379 if (is.null(lop$thiskeyset)) { 380 errex.AFNI("Failed to create thiskeyset from labels"); 381 } 382 } 383 384 if ( !is.null(lop$groupedlabeltablefile) && 385 !file_test('-f', lop$groupedlabeltablefile)) { 386 errex.AFNI(sprintf('labeltable file "%s" does not exist', 387 lop$groupedlabeltablefile)); 388 } 389 390 if (!is.null(lop$Tset) && length(lop$Tset) > 1) { 391 sv <- lop$Tset[[1]]$subj; 392 for (kk in 2:length(lop$Tset)) { 393 sv <- c(sv, lop$Tset[[kk]]$subj) 394 } 395 sv <- unique(sv) 396 if (length(sv) != length(lop$Tset)) { 397 398 errex.AFNI(paste("Have", length(sv), "unique test subject labels (", 399 paste(sv, collapse=','), 400 " ) but", length(lop$Tset)," '-test/-class' options" ) ) 401 402 } 403 } 404 if (!is.null(lop$tset) && length(lop$tset) > 1) { 405 sv <- lop$tset[[1]]$subj; 406 for (kk in 2:length(lop$tset)) { 407 sv <- c(sv, lop$tset[[kk]]$subj) 408 } 409 sv <- unique(sv) 410 if (length(sv) != length(lop$tset)) { 411 412 errex.AFNI(paste("Have", length(sv), "unique test subject labels (", 413 paste(sv, collapse=','), 414 " ) but", length(lop$tset)," '-train' options" ) ) 415 416 } 417 } 418 419 if (lop$method == 'svm' || lop$method == 'SVM') { 420 lop$method <- 'svm' 421 if (is.null(lop$samples_frac)) lop$samples_frac <- 0.03 422 if (is.null(lop$max_samples)) lop$max_samples <- 100 423 if (is.null(lop$min_samples)) lop$min_samples <- 100 424 } else if (lop$method == 'prob' || lop$method == 'PROB') { 425 lop$method <- 'prob' 426 if (is.null(lop$samples_frac)) lop$samples_frac <- 0.03 427 if (is.null(lop$max_samples)) lop$max_samples <- 100 428 if (is.null(lop$min_samples)) lop$min_samples <- 100 429 } 430 return(lop) 431}# end of read.SigsClassify.opts.batch 432 433#Change options list to SigsClassify variable list 434process.SigsClassify.opts <- function (lop, verb = 0) { 435 return(lop) 436} 437 438 439prob.Train.SigsClassify <- function (lvols, samples_frac=NULL, 440 max_samples = Inf, min_samples = 0, 441 thiskeyset=NULL, ltfile=NULL, 442 no_X11=FALSE, prefix="prob.Train.SigsClassify", 443 verb = 1) { 444 if (is.null(lsamp <- get.SampleSet(lvols, samples_frac = samples_frac, 445 max_samples = max_samples, 446 min_samples = min_samples, 447 thiskeyset=thiskeyset, ltfile=ltfile, prefix=prefix, 448 verb = verb, no_X11=no_X11))) { 449 450 errex("Failed to get sample set"); 451 } 452 453 #Copy all fields from lsamp into lprob 454 lprob <- list(sigs = NULL, labs = NULL, Nlabsamp=NULL, 455 model=NULL, labeltable=NULL) 456 for (i in 1:length(lsamp)) { 457 nn <- names(lsamp[i]) 458 lprob[nn] <- lsamp[i] 459 } 460 rm(lsamp) 461 462 if (verb) note.AFNI(sprintf( 463 "Prob model building %d samples of %d classes, and %d features", 464 dim(lprob$sigs)[1], length(lprob$Nlabsamp), dim(lprob$sigs)[2] ), 465 paste(' Label:',names(lprob$Nlabsamp), 466 ', Nsamples',lprob$Nlabsamp, '\n', 467 collapse=''), 468 tic=1) 469 470 471 #foreach feature, and foreach class, fit the distribution 472 lprob$pfci <- matrix( 473 ncol=length(lprob$class.set), 474 nrow=length(lprob$feat.labels)) 475 colnames(lprob$pfci) <- lprob$class.set 476 rownames(lprob$pfci) <- lprob$feat.labels 477 lprob$pfc <- list() 478 for (feat in lprob$feat.labels) { 479 for (cls in lprob$class.set) { 480 note.AFNI(paste("Fitting ", feat, cls)) 481 #You must apply proper scaling (scl=1), or else 482 #fit fails in many instances 483 smp <- get.train.samples(lprob, feat=feat, cls=cls, scl = 1, 484 verb=max(0,verb-1)) 485 disp.train.samples(smp); 486 #prompt.AFNI() 487 kk<-list(fd=fit.pdist(smp)) 488 lprob$pfc <- c(lprob$pfc, kk) 489 lprob$pfci[feat,cls] <- length(lprob$pfc) 490 } 491 } 492 493 494 return(lprob); 495 496} 497 498svm.Train.SigsClassify <- function (lvols, samples_frac=NULL, 499 max_samples = Inf, min_samples = 0, 500 thiskeyset=NULL, ltfile=NULL, 501 no_X11=FALSE, prefix="svm.Train.SigsClassify", 502 svmscale=TRUE, 503 verb = 1, no_tune = FALSE, doprob=TRUE) { 504 if (is.null(lsamp <- get.SampleSet(lvols, samples_frac = samples_frac, 505 max_samples = max_samples, 506 min_samples = min_samples, 507 thiskeyset=thiskeyset, ltfile=ltfile, prefix=prefix, 508 verb = verb))) { 509 510 errex.AFNI("Failed to get sample set"); 511 } 512 #Copy all fields from lsamp into lsvm 513 lsvm <- list(sigs = NULL, labs = NULL, Nlabsamp=NULL, 514 model=NULL, labeltable=NULL) 515 for (i in 1:length(lsamp)) { 516 nn <- names(lsamp[i]) 517 lsvm[nn] <- lsamp[i] 518 } 519 rm(lsamp) 520 521 if (verb) note.AFNI(sprintf( 522 "SVM model building %d samples of %d classes, and %d features", 523 dim(lsvm$sigs)[1], length(lsvm$Nlabsamp), dim(lsvm$sigs)[2] ), 524 paste(' Label:',names(lsvm$Nlabsamp), 525 ', Nsamples',lsvm$Nlabsamp, '\n', 526 collapse=''), 527 tic=1) 528 529 530 #if (verb>1) browser() 531 if (no_tune) { 532 if (verb) note.AFNI(sprintf("NO fine tuning. Probability=%d", doprob)) 533 lsvm$model <- svm(lsvm$sigs, factor(lsvm$labs), 534 scale = svmscale, probability=doprob, 535 type="C-classification") 536 } else { 537 if (verb) note.AFNI(sprintf("With fine tuning. Probability=%d", doprob)) 538 lsvm$modeltuned <- tune(svm, 539 train.x=lsvm$sigs, train.y=factor(lsvm$labs), 540 ranges = list( gamma = c(2^(-6:3), 1/dim(lsvm$sigs)[2]), 541 cost = 2^(-4:4)), 542 scale = svmscale, probability=doprob, 543 type="C-classification" ) 544 lsvm$model <- lsvm$modeltuned$best.model 545 } 546 547 if (verb) note.AFNI(sprintf("SVM modeling completed"), tic=1) 548 549 return(lsvm) 550} 551 552prob.Test.SigsClassify <- function (lvol, lprob, i, ltfile, onames, 553 mixfrac, verb=1) { 554 if (mixfrac=='uniform') { 555 mixfrac <- rep(1,length(lprob$class.set))/length(lprob$class.set) 556 names(mixfrac) <- lprob$class.set 557 } else { 558 err.AFNI(paste("Not ready for mixfrac other than uniform.\nHave ", 559 mixfrac)) 560 return(NULL) 561 } 562 563 #Now do the prob. testing 564 note.AFNI(sprintf("About to PROB test on %s", lvol[[1]]$signame)); 565 #if (verb>1) browser() 566 ivox <- dset.nonzero.indices(lvol[[1]]$sig) 567 568 if (0) { #WAY too slow 569 #Compute probabilities 570 dd <- c(lvol[[1]]$sig$dim[1:3],length(lprob$class.set)) 571 brk <- array(0,dd) 572 dim(brk) <- c(prod(dd[1:3]),dd[4]) 573 uid <- newid.AFNI() 574 uid <- sprintf('subj%d',i) 575 for (v in 1:dd[4]) { 576 note.AFNI(paste("SB ",v)); 577 brk[ivox,v] <- p.cv_GIV_A(lvol[[1]]$sig$brk[ivox,], 578 lprob$feat.labels, cls = lprob$class.set[v], 579 lprob, mixfrac, uid = uid, verb=0) 580 } 581 note.AFNI("Writing probs"); 582 write.AFNI(test.vol(lvol,onames,'p','prob'), brk, 583 label=paste('p', lprob$class.set, sep='.'), 584 defhead=lvol[[1]]$sig$header, 585 verb = max(0,verb-1), scale=TRUE) 586 } else { 587 nn <- parse.AFNI.name(lvol[[1]]$signame) 588 pp <- parse.AFNI.name(test.vol(lvol,onames,'p','prob')) 589 cc <- parse.AFNI.name(test.vol(lvol,onames,'c','prob')) 590 if (exists.AFNI.name(test.vol(lvol,onames,'p','prob'))) { 591 note.AFNI(sprintf("Volume %s exists. Reusing...", 592 test.vol(lvol,onames,'p','prob'))); 593 } else { 594 anatopt <- '' 595 if (!is.null(lvol[[1]]$anatname)) 596 anatopt <- sprintf('-anat %s', lvol[[1]]$anatname) 597 cm <- paste('3dGenPriors -sig ', lvol[[1]]$signame, 598 '-uid ', sprintf('%s%s',lop$uidp, nn$prefix), 599 '-do pc', '-tdist ', sprintf("%s.niml.td",lop$trainsuffix), 600 anatopt, 601 '-labeltable ', ltfile, 602 '-pprefix ', pp$pprefix, 603 '-cprefix ', cc$pprefix) 604 sys.AFNI(cm, echo=TRUE); 605 } 606 #now loadit back in 607 if (is.null(pset <- 608 dset.3DBRKarrayto1D(read.AFNI(test.vol(lvol,onames,'p','prob'))))) { 609 errex.AFNI("Failed to read output"); 610 } 611 } 612 if (0) { #Now done in GenPriors 613 #Do the classes 614 note.AFNI("Getting imax prob"); 615 ddm <- c(lvol[[1]]$sig$dim[1:3],1) 616 brkm <- array(0,ddm) 617 dim(brkm) <- c(prod(ddm[1:3]),ddm[4]) 618 for (i in ivox) { 619 brkm[i] <- lprob$classkey.set[which.max(pset$brk[i,])] 620 } 621 dim(brkm) <- ddm 622 note.AFNI("Writing max prob"); 623 write.AFNI( test.vol(lvol,onames,'c','prob'), brkm, 624 label=c('maxprob_labels'), defhead=lvol[[1]]$sig$header, 625 verb = max(0,verb-1), scale=FALSE) 626 if (!is.null(ltfile)) { 627 scom <- sprintf('3drefit -labeltable %s %s', 628 ltfile, test.vol(lvol,onames,'c','prob')); 629 system(scom) 630 } 631 } 632 return() 633} 634 635svm.Test.SigsClassify <- function (lvol,lsvm,i, doprob, ltfile, onames, verb=1) { 636 #Now do the SVM testing 637 note.AFNI(sprintf("About to SVM test on %s, probability = %d", 638 lvol[[1]]$signame, doprob)); 639 #if (verb>1) browser() 640 ivox <- dset.nonzero.indices(lvol[[1]]$sig) 641 classy<-predict(lsvm$model, lvol[[1]]$sig$brk[ivox,], 642 probability=doprob, 643 decision.values = TRUE) 644 645 646 #Make a new volume data 647 dd <- c(lvol[[1]]$sig$dim[1:3],1) 648 brk <- array(0,dd) 649 #Change brk to 1D 650 dim(brk) <- c(prod(dd[1:3]),dd[4]) 651 brk[ivox,1] <- 652 as.integer(levels(classy))[classy] 653 #Fastest way to get factor levels back as integer labels 654 #http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg08378.html 655 #brk back to 3D 656 dim(brk) <- dd 657 658 #write out dset 659 write.AFNI( test.vol(lvol,onames,'c','svm'), brk, 660 label=c('test_labels'), defhead=lvol[[1]]$sig$header, 661 verb = max(0,verb-1), scale=FALSE) 662 663 if (!is.null(ltfile)) { 664 scom <- sprintf('3drefit -labeltable %s %s', 665 ltfile, test.vol(lvol,onames,'c','svm')); 666 system(scom) 667 } 668 669 dv <- as.matrix(attr(classy, "decision.values")) 670 dd <- c(lvol[[1]]$sig$dim[1:3],dim(dv)[2]) 671 brk <- array(0,dd) 672 dim(brk) <- c(prod(dd[1:3]),dd[4]) 673 brk[ivox,] <- dv 674 dv<-NULL 675 dim(brk) <- dd 676 write.AFNI(test.vol(lvol,onames,'dv','svm'), brk, 677 label=paste('tdecis_val', seq(1,dd[4]), sep='.'), 678 defhead=lvol[[1]]$sig$header, 679 verb = max(0,verb-1), scale=TRUE) 680 681 if (doprob) { 682 dv <- as.matrix(attr(classy, "probabilities")) 683 dd <- c(lvol[[1]]$sig$dim[1:3],dim(dv)[2]) 684 brk <- array(0,dd) 685 dim(brk) <- c(prod(dd[1:3]),dd[4]) 686 brk[ivox,] <- dv 687 dv<-NULL 688 dim(brk) <- dd 689 if (dd[4] != length(lsvm$class.set)) { 690 err.AFNI(paste ("have", dd[4], "subbricks, but ", 691 length(lsvm$labeltable$labels), "training labels\n", 692 "Something is fishy here")); 693 return(NULL); 694 } 695 write.AFNI(test.vol(lvol,onames,'p','svm'), brk, 696 label=paste('p', lsvm$class.set, sep='.'), 697 defhead=lvol[[1]]$sig$header, 698 verb = max(0,verb-1), scale=TRUE) 699 } 700 701 #Need space so delete a little 702 if (verb > 1) { 703 llvv <- lvol 704 llvv <- c(llvv, classy=classy, ivox=ivox) 705 fout <- sprintf('%s.lvols%d.Rdat',test.vol(lvol,onames,'c','svm'),i) 706 note.AFNI(sprintf("Saving svm result in %s", fout)); 707 save(llvv, file=fout) 708 llvv <- NULL 709 } 710 711 return() 712} 713 714#Standard names of testing output volumes 715test.vol <- function (lvol, onames, tp, mode) { 716 volsuffix <- onames$volsuffix 717 opath <- onames$opath 718 opref <- onames$opref 719 720 if (is.null(volsuffix)) volsuffix <- '' 721 else { 722 #make sure it has a . 723 if (length(grep ('^\\.', volsuffix)) == 0) { 724 volsuffix <- sprintf('.%s', volsuffix); 725 } 726 } 727 if (!is.null(opath)) { 728 if (!file_test('-d', opath)) { 729 err.AFNI(sprintf('Output path %s does not exist', opath)); 730 return(NULL); 731 } 732 } 733 #prep names 734 an <- parse.AFNI.name(lvol[[1]]$signame) 735 if (is.null(opath)) ppo <- an$path 736 else ppo <- opath 737 if (is.null(opref)) { 738 ppp <- sprintf('%s.%s%s', lvol[[1]]$subj, an$prefix, volsuffix) 739 } else { 740 ppp <- sprintf('%s.%s%s', lvol[[1]]$subj, opref, volsuffix) 741 } 742 if (mode=='svm') { mm <- 'S' } 743 else if (mode=='prob') { mm <- 'P' } 744 else return(NULL) 745 746 if (tp == 'c') return(sprintf('%s/%s.%s.cls%s' , ppo, ppp, mm, an$view)) 747 else if (tp == 'p') return(sprintf('%s/%s.%s.clsp%s' , ppo, ppp, mm, an$view)) 748 else if (tp == 'dv') return(sprintf('%s/%s.%s.clsdv%s', 749 ppo, ppp, mm, an$view)) 750 else if (tp == 'dc') return(sprintf('%s/%s.%s.DiceCoef', ppo, ppp, mm)) 751 else if (tp == 'gc') return(sprintf('%s/%s.%s.gcls%s' , 752 ppo, ppp, mm, an$view)) 753 else if (tp == 'gp') return(sprintf('%s/%s.%s.gclsp%s' , 754 ppo, ppp, mm, an$view)) 755 else if (tp == 'gdv') return(sprintf('%s/%s.%s.gclsdv%s', 756 ppo, ppp, mm, an$view)) 757 else if (tp == 'gdc') return(sprintf('%s/%s.%s.gDiceCoef', ppo, ppp, mm)) 758 else if (tp == 'gld') return(sprintf('%s/%s.%s.cls%s', 759 ppo, lvol[[1]]$subj, 'gold',an$view)) 760 else if (tp == 'ggld') return(sprintf('%s/%s.%s.gcls%s', 761 ppo, lvol[[1]]$subj, 'gold',an$view)) 762 else if (tp == 'a') return(sprintf('%s/%s.%s%s', 763 ppo, lvol[[1]]$subj, 'anat',an$view)) 764 else return(NULL) 765 766} 767 768Test.SigsClassify <- function (lvols, ltrain=NULL, verb = 1, 769 ltfile=NULL, overwrite=FALSE, 770 onames = list(volsuffix=NULL, 771 opath = NULL, opref = NULL), 772 ltgroupedfile = NULL, doprob=TRUE, 773 method='svm', mixfrac='uniform') { 774 775 if (method=='svm') { 776 if (is.null(ltrain) || is.null(ltrain$model)) { 777 err.AFNI("Bad input") 778 return(NULL); 779 } 780 } else if (method == 'prob') { 781 if (is.null(ltrain) ) { 782 err.AFNI("Bad input") 783 return(NULL); 784 } 785 } else return(NULL); 786 787 for (i in 1:length(lvols)) { 788 789 if (!overwrite && 790 exists.AFNI.name(test.vol(lvols[i],onames,'c',method))) { 791 note.AFNI(sprintf("Volume %s exists. Skipping...", 792 test.vol(lvols[i],onames,'c',method))); 793 } else { 794 if (verb) note.AFNI(sprintf("About to process signatures %s", 795 lvols[i][[1]]$signame)); 796 #Load signatures 797 if (is.null(lvols[i][[1]]$sig)) { 798 if (verb > 1) note.AFNI(sprintf("Loading %s", 799 lvols[i][[1]]$signame),tic=1) 800 lvols[i][[1]]$sig <- 801 dset.3DBRKarrayto1D(read.AFNI(lvols[i][[1]]$signame)) 802 } else { 803 if (verb > 1) note.AFNI(sprintf("signatures already loaded")) 804 lvols[i][[1]]$sig <- dset.3DBRKarrayto1D(lvols[i][[1]]$sig) 805 } 806 if (is.null(lvols[i][[1]]$sig)) { 807 err.AFNI("NULL sig"); 808 return(NULL); 809 } 810 811 812 lvols[i][[1]] <- c(lvols[i][[1]], 813 list(classout=test.vol(lvols[i],onames,'c',method))) 814 if (method == 'svm') { 815 svm.Test.SigsClassify(lvols[i], ltrain, i, doprob, ltfile, 816 onames, verb) 817 } else { 818 prob.Test.SigsClassify(lvols[i], ltrain, i, ltfile, onames, 819 mixfrac, verb) 820 } 821 #delete the big stuff. everything's been pretty much saved 822 #This can grow quite large with probabilities 823 if (verb) memory.hogs(msg='Pre NULLing') 824 lvols[i][[1]]$sig$brk <- NULL; 825 gc(); 826 if (verb) memory.hogs(msg='Post NULLing') 827 828 } 829 830 #Put a copy of the anat in place 831 if (!is.null(lvols[i][[1]]$anatname)) { 832 if (!copy.AFNI.dset(lvols[i][[1]]$anatname, 833 test.vol(lvols[i],onames,'a',method), overwrite = TRUE)) { 834 err.AFNI(sprintf( 835 "Failed to create copy of anatomical volume %s", 836 test.vol(lvols[i],onames,'a',method))) 837 return(NULL); 838 } 839 } 840 841 #Create a merged version if necessary 842 if (!is.null(ltgroupedfile) && 843 ( overwrite || 844 !exists.AFNI.name(test.vol(lvols[i],onames,'gc',method)) || 845 !exists.AFNI.name(test.vol(lvols[i],onames,'gp',method))) ) { 846 if (exists.AFNI.name(test.vol(lvols[i],onames,'gc',method))) 847 system(sprintf('\\rm -f %s.HEAD %s.BRIK*', 848 test.vol(lvols[i],onames,'gc',method), 849 test.vol(lvols[i],onames,'gc',method))) 850 if (doprob && 851 exists.AFNI.name(test.vol(lvols[i],onames,'gp',method))) 852 system(sprintf('\\rm -f %s.HEAD %s.BRIK*', 853 test.vol(lvols[i],onames,'gp',method), 854 test.vol(lvols[i],onames,'gp',method))) 855 856 if (doprob) { 857 if (group.vollabels( 858 labelvol = test.vol(lvols[i],onames,'c',method), 859 pvol=test.vol(lvols[i],onames,'p',method), 860 grptable=ltgroupedfile, 861 grplabelpref=test.vol(lvols[i],onames,'gc',method), 862 grpppref=test.vol(lvols[i],onames,'gp',method)) == 0) { 863 err.AFNI("Failed to group labels"); 864 return(NULL); 865 } 866 } else { 867 if (group.vollabels( 868 labelvol = test.vol(lvols[i],onames,'c',method), 869 pvol=NULL, 870 grptable=ltgroupedfile, 871 grplabelpref=test.vol(lvols[i],onames,'gc',method), 872 grpppref=NULL) == 0) { 873 err.AFNI("Failed to group labels"); 874 return(NULL); 875 } 876 } 877 878 } 879 #Do some dice action 880 if (!is.null(lvols[i][[1]]$labname)) { 881 if (!copy.AFNI.dset(lvols[i][[1]]$labname, 882 test.vol(lvols[i],onames,'gld',method), 883 overwrite = TRUE)) { 884 err.AFNI(sprintf( 885 "Failed to create copy of reference label volume %s", 886 test.vol(lvols[i],onames,'gld',method))) 887 return(NULL); 888 } 889 lvols[i][[1]]$dice <- 890 dice.vollabels( 891 basevol=test.vol(lvols[i],onames,'gld',method), 892 labeltable=ltfile, 893 prefix=test.vol(lvols[i],onames,'dc',method), 894 labelvol=test.vol(lvols[i],onames,'c',method), 895 maskby ='LABEL') 896 if (!is.null(ltgroupedfile)) { 897 #first group the base 898 if (group.vollabels( 899 labelvol = test.vol(lvols[i],onames,'gld',method), 900 pvol=NULL, 901 grptable=ltgroupedfile, 902 grplabelpref=test.vol(lvols[i],onames,'ggld',method), 903 grpppref=NULL) == 0) { 904 err.AFNI("Failed to group base labels"); 905 return(NULL); 906 } 907 #Now do the dice 908 lvols[i][[1]]$gdice <- dice.vollabels( 909 test.vol(lvols[i],onames,'ggld',method), 910 labeltable=ltgroupedfile, 911 prefix=test.vol(lvols[i],onames,'gdc',method), 912 labelvol=test.vol(lvols[i],onames,'gc',method), 913 maskby ='LABEL') 914 } 915 } 916 917 } 918 919 return(lvols) 920} 921 922 923################################################################################# 924########################## Begin SigsClassify main ################################### 925################################################################################# 926 927 928 if (!exists('.DBG_args')) { 929 args = (commandArgs(TRUE)) 930 rfile <- first.in.path(sprintf('%s.R',ExecName)) 931 farg <- sprintf("%s.dbg.AFNI.args", ExecName) 932 # save only on -dbg_args 28 Apr 2016 [rickr] 933 if ( '-dbg_args' %in% args ) { 934 save(args, rfile, file=farg, ascii = TRUE) 935 note.AFNI(paste("Saved command line arguments in ", farg)) 936 } 937 } else { 938 note.AFNI("Using .DBG_args resident in workspace"); 939 args <- .DBG_args 940 } 941 if (!length(args)) { 942 BATCH_MODE <<- 0 943 lop <- init.SigsClassify.lop() 944 lop$verb=1 945 note.AFNI("I am not going to have a hard coded list of options here. 946 For testing, run the command line version once. Then launch R in that 947 directory and do the following: 948 source('~/AFNI/src/R_scripts/AFNIio.R'); 949 load.debug.AFNI.args() 950 source('~/AFNI/src/R_scripts/3dSignatures.R');"); 951 } else { 952 if (!exists('.DBG_args')) { 953 BATCH_MODE <<- 1 954 } else { 955 BATCH_MODE <<- 0 956 } 957 if (is.null(lop <- read.SigsClassify.opts.batch(args, verb = 0))) { 958 stop('Error parsing input'); 959 } 960 #str(lop); 961 if (is.null(lop <- process.SigsClassify.opts(lop, verb = lop$verb))) { 962 stop('Error processing input'); 963 } 964 } 965 if (lop$verb) { 966 str(lop); 967 } 968 969 if (is.null(lop$method) || lop$method != 'prob') { 970 note.AFNI("SVM approach"); 971 if (!is.null(lop$load_train)) { 972 if (!file.exists(lop$load_train)) { 973 974 if (file.exists(kk <- sprintf('%s.Rdat',lop$load_train))) { 975 lop$load_train <- kk 976 } else { 977 errex.AFNI(c(lop$load_train, " does not exist.")); 978 } 979 } 980 981 note.AFNI(sprintf("Loading training set from pre-existing %s", 982 lop$load_train )) 983 if (exists('train')) { 984 rm('train') 985 } 986 #Cannot override existing lop !, some .Rdat files 987 #contain lop <- stupid thing to do 988 lop_buf <- lop; rm('lop'); 989 load(lop_buf$load_train); 990 if (!exists('train') || is.null(train)) { 991 warn.AFNI("Variable 'train' not found or null!") 992 train <- NULL 993 } 994 if (exists('lop')) {#part of the bug fix above 995 lop_train <- lop; 996 rm('lop'); 997 } 998 lop <- lop_buf; rm('lop_buf'); 999 } else { 1000 if (is.null(lop$trainsuffix)) { 1001 errex.AFNI("Have no trainsuffix"); 1002 } 1003 trainname <- sprintf("%s.Rdat",lop$trainsuffix); 1004 train <- NULL 1005 if (lop$reuse_train) { 1006 if (file.exists(trainname)) { 1007 note.AFNI(sprintf("Loading training set from pre-existing %s", 1008 trainname)) 1009 rm('train') 1010 lop_buf <- lop; rm('lop'); 1011 load(trainname); 1012 if (!exists('train') || is.null(train)) { 1013 warn.AFNI("Variable 'train' not found or null!") 1014 train <- NULL 1015 } 1016 if (exists('lop')) {#part of the bug fix above 1017 lop_train <- lop; 1018 rm('lop'); 1019 } 1020 lop <- lop_buf; rm('lop_buf'); 1021 } 1022 } 1023 } 1024 if (is.null(train)) { 1025 if (length(lop$tset) < 1) 1026 errex.AFNI("Have no training data or training datasets") 1027 note.AFNI("Creating new training object"); 1028 train <- svm.Train.SigsClassify(lop$tset, samples_frac=lop$samples_frac, 1029 max_samples = lop$max_samples, 1030 min_samples = lop$min_samples, 1031 thiskeyset=lop$thiskeyset, 1032 ltfile = lop$labeltablefile, 1033 prefix=lop$trainsuffix, 1034 verb=lop$verb, no_X11=lop$no_X11, 1035 svmscale=lop$svmcale, no_tune = lop$no_tune, 1036 doprob = lop$doprob) 1037 if (is.null(train)) errex.AFNI("Failed to get training set"); 1038 lop_train <- lop 1039 save(train, lop_train, 1040 file=sprintf("%s.Rdat",lop$trainsuffix), ascii = TRUE) 1041 } else { 1042 note.AFNI("Reusing/Or reloading SVM train object"); 1043 if (is.null(train$labeltable)) { 1044 uk <- unique(train$labs) 1045 warn.AFNI(paste("Old format training struct.\n", 1046 "Assuming ", lop$labeltablefile, 1047 "is the proper label table for label keys [", 1048 paste(uk,collapse=' '), "]", 1049 collapse='')); 1050 1051 train$labeltable <- read.AFNI.labeltable(ltfile=lop$labeltablefile) 1052 if (is.null(train$labeltable)) { 1053 err.AFNI("Failed to build test label table list"); 1054 return(NULL) 1055 } 1056 } 1057 1058 } 1059 if (is.null(train)) { 1060 errex.AFNI("End of the line, no train here"); 1061 } 1062 #browser() 1063 test <- Test.SigsClassify(lop$Tset, ltrain=train, 1064 verb=lop$verb, ltfile = lop$labeltablefile, 1065 onames = list(volsuffix = lop$volsuffix, 1066 opath = lop$odir, opref = lop$opref), 1067 ltgroupedfile = lop$groupedlabeltablefile, 1068 doprob = lop$doprob, method = lop$method) 1069 lop_test <- lop 1070 save(test, lop_test, file=sprintf("%s.Rdat",lop$testsuffix), ascii = TRUE) 1071 } else { 1072 #Prob. based approach 1073 if (!is.null(lop$load_train)) { 1074 if (!file.exists(lop$load_train)) { 1075 if (file.exists(kk <- sprintf('%s.Rdat',lop$load_train))) { 1076 lop$load_train <- kk 1077 } else { 1078 errex.AFNI(paste(lop$load_train, " does not exist.")); 1079 } 1080 } 1081 1082 note.AFNI(sprintf("Loading training set from pre-existing %s", 1083 lop$load_train )) 1084 if (exists('train')) rm('train') 1085 1086 load(lop$load_train); 1087 if (!exists('train') || is.null(train)) { 1088 warn.AFNI("Variable 'train' not found or null!") 1089 train <- NULL 1090 } 1091 } else { 1092 if (is.null(lop$trainsuffix)) { 1093 errex.AFNI("Have no trainsuffix"); 1094 } 1095 trainname <- sprintf("%s.Rdat",lop$trainsuffix); 1096 train <- NULL 1097 if (lop$reuse_train) { 1098 if (file.exists(trainname)) { 1099 note.AFNI(sprintf("Loading training set from pre-existing %s", 1100 trainname)) 1101 rm('train') 1102 load(trainname); 1103 if (!exists('train') || is.null(train)) { 1104 warn.AFNI("Variable 'train' not found or null!") 1105 train <- NULL 1106 } 1107 } 1108 } 1109 } 1110 if (is.null(train)) { 1111 if (length(lop$tset) < 1) 1112 errex.AFNI("Have no training data or training datasets") 1113 note.AFNI("Creating new training object"); 1114 train <- prob.Train.SigsClassify( 1115 lop$tset, samples_frac=lop$samples_frac, 1116 max_samples = lop$max_samples, 1117 min_samples = lop$min_samples, 1118 prefix=lop$trainsuffix, 1119 thiskeyset=lop$thiskeyset, 1120 ltfile = lop$labeltablefile, 1121 verb=lop$verb, no_X11=lop$no_X11) 1122 if (is.null(train)) errex.AFNI("Failed to get training set"); 1123 lop_train <- lop 1124 save(train, lop_train, 1125 file=sprintf("%s.Rdat",lop$trainsuffix), ascii = TRUE) 1126 NI.write.train.dists(train, fname=sprintf("%s.niml.td",lop$trainsuffix)) 1127 } else { 1128 note.AFNI("Reusing/Or reloading PROB train object"); 1129 } 1130 if (is.null(train)) { 1131 errex.AFNI("End of the line, no train here"); 1132 } 1133 #browser() 1134 test <- Test.SigsClassify(lop$Tset, ltrain=train, 1135 verb=lop$verb, ltfile = lop$labeltablefile, 1136 onames= list(volsuffix = lop$volsuffix, 1137 opath = lop$odir, opref = lop$opref), 1138 ltgroupedfile = lop$groupedlabeltablefile, 1139 method=lop$method, mixfrac='uniform') 1140 lop_test <- lop 1141 save(test, lop_test, file=sprintf("%s.Rdat",lop$testsuffix), ascii = TRUE) 1142 1143 } 1144 note.AFNI("All done."); 1145