1library(lattice) 2library(fitdistrplus) 3 4#------------------------------------------------------------------ 5# Functions to support 3dSignatures.R 6#------------------------------------------------------------------ 7groupofKey.labeltable <- function (key, labeltable) { 8 if (is.character(labeltable)) 9 labeltable <- read.AFNI.labeltable(labeltable) 10 gok <- vector() 11 for (k in key) { 12 g <- which (labeltable$keys == k) 13 if (length(g)) gok <- c(gok,labeltable$groups[g]) 14 else { 15 err.AFNI("Failed") 16 return(NULL) 17 } 18 } 19 return(gok) 20} 21 22keylabel.labeltable <- function (key=NULL, ltfile=NULL) { 23 if (0) { #OLD 24 com <- sprintf( 25 '@MakeLabelTable -quiet_death -labeltable %s -klabel %d', 26 ltfile, key) 27 lab <- system(com, ignore.stderr = TRUE, intern=TRUE) 28 if (length(lab) == 0) { 29 err.AFNI(paste("Failed to get key label from labeltable.\n", 30 "Command ",com,"Failed")); 31 return(NULL); 32 } 33 } else { 34 if (is.character(ltfile)) 35 ltfile <- read.AFNI.labeltable(ltfile) 36 lab <- vector() 37 for (kk in key) { 38 lab <-c(lab, ltfile$labels[which(ltfile$keys==kk)]) 39 } 40 } 41 return(lab); 42} 43 44labelkey.labeltable <- function (label=NULL, ltfile=NULL) { 45 if (0) { #old 46 kk <- vector() 47 for (ll in label) { 48 com <- sprintf( 49 '@MakeLabelTable -word_label_match -quiet_death -labeltable %s -lkeys %s', 50 ltfile, ll) 51 lab <- system(com, ignore.stderr = TRUE, intern=TRUE) 52 if (length(lab) == 0) { 53 err.AFNI(paste("Failed to get key of label", 54 label," from labeltable",ltfile, ".\n", 55 "Command ",com,"Failed")); 56 return(NULL); 57 } 58 kk <- c(kk, as.numeric(lab)) 59 60 if (is.na(kk[length(kk)])) { 61 err.AFNI(paste("Failed to get key of label", 62 label," from labeltable.\n", 63 "Command ",com,"Failed")); 64 return(NULL); 65 } 66 } 67 } else { 68 if (is.character(ltfile)) 69 ltfile <- read.AFNI.labeltable(ltfile) 70 kk <- vector() 71 for (ll in label) { 72 kk <- c(kk, ltfile$keys[which(ltfile$labels==ll)]) 73 } 74 75 } 76 return(kk) 77} 78 79build.labeltable <- function(labeltable=NULL, keys=NULL, 80 grp.label = c('CSF','GM','WM','InSk','OutSk','Out')) { 81 #OBSOLETE, use read.AFNI.labeltable instead 82 if (is.null(keys)) { 83 com <- sprintf( 84 '@MakeLabelTable -quiet_death -labeltable %s -all_keys', labeltable) 85 ss <- sys.AFNI(com) 86 if (ss$stat == 1) { 87 err.AFNI(paste("Failed to get keys from labeltable",ltfile, ".\n", 88 "Command ",com,"Failed")); 89 return(NULL); 90 } 91 keys <- as.numeric(ss$out) 92 } 93 lt <- list(labels=NULL, keys=keys, grp.label=grp.label, groups=NULL) 94 for (key in keys) { 95 if (is.null(lab <- keylabel.labeltable(key=key, ltfile=labeltable))) { 96 return(NULL); 97 } 98 99 lt$labels <- c(lt$labels, lab) 100 if (!is.null(grp.label)) { 101 nn <- 0 102 ig <- 0 103 for (i in 1:length(grp.label)) { 104 if (length(grep(grp.label[i],lab))) { 105 if (ig == 0) { 106 nn <- 1 107 ig <- i 108 } else { 109 if ( length(grp.label[ig]) != length(lab) ) { 110 if (length(grp.label[i]) == length(lab) ) { 111 # A better match 112 nn <- 1 113 ig <- i 114 } else { 115 #Another partial 116 nn <- nn + 1 117 ig <- i 118 } 119 } 120 } 121 } 122 123 } 124 if (nn != 1) { 125 warn.AFNI(paste( 126 "Failed to find 1 group for ",lab, ". Found ", 127 nn ,"candidates. Setting grp to", 128 length(grp.label)+1)) 129 ig <- length(grp.label)+1 130 } 131 lt$groups <- c(lt$groups, ig) 132 133 } 134 } 135 if (!is.null(grp.label)) { 136 if (max(lt$groups) > length(grp.label)) { 137 lt$grp.label <- c(grp.label,"Unknown") 138 } else { 139 lt$grp.label <- grp.label 140 } 141 } 142 return(lt) 143} 144 145labels.labeltable <- function(ltfile=NULL) { 146 if (0) { #OLD APPROACH 147 com <- paste('@MakeLabelTable -labeltable ', ltfile, 148 '-all_labels', 149 collapse = ' '); 150 labs <- system(com, ignore.stderr = TRUE, intern=TRUE) 151 if (length(labs) == 0) { 152 err.AFNI(paste("Failed to get key label from labeltable.\n", 153 "Command ",com,"Failed")); 154 return(NULL); 155 } 156 return(labs) 157 } else { 158 if (is.character(ltfile)) 159 ltfile <- read.AFNI.labeltable(ltfile) 160 return(ltfile$labels) 161 } 162} 163 164 165dice.vollabels <- function(basevol=NULL, 166 labeltable=NULL, 167 prefix=NULL, labelvol=NULL, 168 diceout=NULL, maskby = 'BASE') { 169 if (is.null(diceout)) { 170 diceout <- sprintf('%s.1D', prefix) 171 } 172 if (is.null(labelvol)) { 173 err.AFNI("Need labelvol") 174 return(NULL) 175 } 176 if (is.null(maskby) || maskby == '') { 177 mm <- '' 178 } else if (maskby == 'BASE') { 179 mm <- '-mask_by_base' 180 } else if (maskby == 'LABEL') { 181 mm <- '-mask_by_dset_vals' 182 } else { 183 err.AFNI(paste("Bad value ", maskby)) 184 return(NULL) 185 } 186 com <- paste('@DiceMetric -base ', basevol, 187 '-save_match -save_diff ', mm, 188 '-prefix ', prefix, 189 '-labeltable ',labeltable , 190 '-dsets ', labelvol, 191 '> ', diceout, 192 collapse = ' ') 193 194 sys.AFNI(com, echo = TRUE) 195 #system(sprintf("echo '#Command: %s' >> %s", com, diceout)); 196 dd <- read.AFNI.matrix(diceout) 197 ll <- list(command=com, coef=dd[,2]) 198 names(ll$coef) <- labels.labeltable(labeltable) 199 200 return(ll) 201} 202 203group.vollabels <- function(labelvol=NULL, pvol=NULL, 204 grptable=NULL, 205 grplabelpref=NULL, grpppref=NULL) { 206 if (!is.null(pvol)) { 207 com <- paste('@RegroupLabels -grouped_labeltable ', 208 grptable, '-labdset ',labelvol , 209 '-grpd_labprefix ', grplabelpref, 210 '-pdset', pvol, 211 '-grpd_pprefix', grpppref, 212 collapse = '') 213 } else { 214 com <- paste('@RegroupLabels -grouped_labeltable ', 215 grptable, '-labdset ',labelvol , 216 '-grpd_labprefix ', grplabelpref, 217 collapse = '') 218 } 219 note.AFNI(com) 220 mm <- system(com, ignore.stderr = FALSE, intern=TRUE) 221 note.AFNI(mm) 222 223 if (is.null(pvol)) return(used.AFNI.prefix(grplabelpref)) 224 else return(used.AFNI.prefix(grplabelpref)*used.AFNI.prefix(grpppref)) 225} 226 227get.SampleSet <- function (lvols, samples_frac=NULL, 228 max_samples = Inf, min_samples = 0, 229 thiskeyset=NULL, ltfile=NULL, prefix=NULL, 230 verb = 1, no_X11=TRUE) { 231 if (length(samples_frac) != 1 || samples_frac > 1 || samples_frac <= 0) { 232 err.AFNI("Bad sample fraction") 233 return(NULL); 234 } 235 if (is.null(ltfile)) { 236 err.AFNI("ltfile is now mandatory") 237 return(NULL) 238 } 239 240 241 lsamp <- list(sigs = NULL, labs = NULL, Nlabsamp=NULL, 242 model=NULL, labeltable=NULL, feat.labels=NULL) 243 for (lvol in lvols) { 244 #Load label thingy 245 if (is.null(lvol$lab)) { 246 if (verb) 247 note.AFNI(sprintf("Loading labels volume %s", lvol$labname),tic=1) 248 lvol$lab <- dset.3DBRKarrayto1D(read.AFNI(lvol$labname)) 249 } else { 250 if (verb) note.AFNI(sprintf("Labels volume already loaded")) 251 lvol$lab <- dset.3DBRKarrayto1D(lvol$lab) 252 } 253 if (is.null(lvol$lab)) { 254 err.AFNI('Failed to get labels'); 255 return(NULL) 256 } 257 if (is.null(lvol$labselecs)) { 258 if (verb > 1) note.AFNI(sprintf("Selecting set"),tic=1) 259 if (is.null(thiskeyset)) { 260 if (is.null(classkey.set)) 261 classkey.set <- setdiff(unique(lvol$lab$brk), 0) 262 } else { 263 classkey.set <- setdiff(thiskeyset, 0) 264 } 265 #name classkey.set by the labels from the labeltable 266 lsamp$labeltable <- read.AFNI.labeltable(ltfile) 267 if (is.null(lsamp$labeltable)) { 268 err.AFNI("Failed to load label table "); 269 return(NULL) 270 } 271 lsamp$classkey.set <- classkey.set 272 lsamp$class.set <- keylabel.labeltable(classkey.set, lsamp$labeltable) 273 274 if (verb > 1) note.AFNI(sprintf("Forming labselecs"),tic=1) 275 labselecs <- list() 276 for (lv in classkey.set) { 277 ilab <- which(lvol$lab$brk == lv) 278 nsmax <- length(ilab); 279 if (!is.null(samples_frac)) { 280 ns <- as.integer(samples_frac*nsmax) 281 if (ns > max_samples) ns <- max_samples 282 else if (ns < min_samples) ns <- min(min_samples, nsmax) 283 ii <- sample(1:nsmax, ns) 284 ilab <- ilab[ii] 285 } 286 if (verb > 1) 287 note.AFNI( 288 sprintf( 289 " Got %d labels (%.2f%%) from max. of %d, for subclass %d", ns, ns/nsmax*100, nsmax, lv),tic=1) 290 labselecs <- 291 c (labselecs, 292 list( list( 293 labval=lv, ilab = ilab) ) ) 294 names(labselecs[length(labselecs)]) <- 295 sprintf('VoxSelection.%d', lv) 296 } 297 } else { 298 errex.AFNI("Not implemented yet"); 299 } 300 #Load signatures 301 if (is.null(lvol$sig)) { 302 if (verb) 303 note.AFNI(sprintf("Loading signatures %s", lvol$signame),tic=1) 304 lvol$sig <- dset.3DBRKarrayto1D(read.AFNI(lvol$signame)) 305 306 } else { 307 lvol$sig <- dset.3DBRKarrayto1D(lvol$sig) 308 } 309 if (is.null(lsamp$feat.labels)) 310 lsamp$feat.labels <- dset.labels(lvol$sig) 311 312 if (is.null(lvol$sig)) { 313 err.AFNI('Failed to get signatures'); 314 return(NULL) 315 } 316 #Now build the set of features for each label 317 if (verb) note.AFNI(sprintf("Building sigs and labs"),tic=1) 318 for (lbsel in labselecs) { 319 if (length(lbsel$ilab)) { 320 lsamp$sigs <- rbind(lsamp$sigs,lvol$sig$brk[lbsel$ilab,]) 321 lsamp$labs <- c(lsamp$labs,rep(lbsel$labval,length(lbsel$ilab))) 322 iul <- which(names(lsamp$Nlabsamp) == as.character(lbsel$labval)) 323 if (length(iul)) { 324 lsamp$Nlabsamp[iul] <- lsamp$Nlabsamp[iul]+length(lbsel$ilab) 325 } else { 326 lsamp$Nlabsamp <- c(lsamp$Nlabsamp, length(lbsel$ilab)) 327 names(lsamp$Nlabsamp)[length(lsamp$Nlabsamp)] <- 328 as.character(lbsel$labval) 329 } 330 } 331 } 332 } 333 334 #Before testing, a safety check 335 mm <- matrix(0, dim(lsamp$sigs)[2], length(lsamp$class.set)) 336 for (l in 1:length(lsamp$class.set)) { 337 for (i in 1:dim(lsamp$sigs)[2]) { 338 ifnd <- which (lsamp$labs==lsamp$class.set[l]) 339 if (length(ifnd)) { 340 mm[i,l] <-median(lsamp$sigs[ifnd,i]) 341 } 342 } 343 } 344 345 plot.1D(dmat = mm, oneplot=1, col.ystack=FALSE, 346 col.grp=groupofKey.labeltable(lsamp$class.set, lsamp$labeltable), 347 grp.label=lsamp$labeltable$grp.label, 348 prefix =sprintf('%s.jpg',prefix), ttl.main=prefix, 349 nodisp = no_X11, save.Rdat=TRUE) 350 351 #And show the sample set 352 note.AFNI("About to show all samples...") 353 disp.all.train.samples( lsamp, nodisp = no_X11, 354 prefix =sprintf('%s.trset.pdf',prefix)) 355 356 #compute the scaling 357 lsamp$feat.scale <- compute.feature.scaling (lsamp) 358 359 return(lsamp) 360} 361 362disp.all.train.samples <- function (lsamp, nodisp = FALSE, prefix='hello.pdf') { 363 364 n_col <- ncol(lsamp$sigs) 365 n_samps <- nrow(lsamp$sigs) 366 367 #Class factor 368 cls <- matrix(ncol=n_col, nrow=n_samps) 369 for (i in 1:n_col) cls[,i] = lsamp$labs 370 cls.f <- factor(cls, levels=lsamp$labeltable$keys, 371 labels=lsamp$labeltable$labels) 372 373 #feature factor 374 feat <- matrix(ncol=n_col, nrow=n_samps) 375 for (i in 1:n_col) feat[,i] = rep(i, n_samps) 376 ll <- gsub('_mm','',lsamp$feat.labels) 377 ll <- gsub('MEDIAN.','m',ll) 378 ll <- gsub('MAD.','M',ll) 379 ll <- gsub('P2skew.','s',ll) 380 feat.f <- factor(feat, levels=seq(1:n_col), 381 labels=ll) 382 383 #setup device 384 P <- list(prefix=prefix, nodisp=nodisp, dev.this=NULL, dev.new=FALSE) 385 P$dev.this <- plot.1D.setupdevice(P) 386 print(xyplot(lsamp$sigs~feat.f|cls.f, scales=list(rot=90))) 387 plot.1D.unsetupdevice(P) 388 389 return(0) 390} 391 392get.train.samples <- function (lsamp, feat=NULL, cls=NULL, scl = 0, verb=0) { 393 #Safety check 394 if (ncol(lsamp$sigs) != length(lsamp$feat.labels)) { 395 err.AFNI(paste("Bad news! You must have ", 396 ncol(lsamp$sigs), "elements in lsamp$feat.labels", 397 "Currently have", length(lsamp$feat.labels))) 398 return(NULL); 399 } 400 401 if (verb) 402 note.AFNI(paste("Fetching samples of feature", feat,", class ", cls)) 403 404 icol <- which(lsamp$feat.labels == feat) 405 if (length(icol) != 1) { 406 err.AFNI(paste("Feature", feat,"not found")) 407 return(NULL); 408 } 409 if (verb > 1) note.AFNI(paste("Column ",icol, "matches feature", feat)) 410 411 if (cls == 'ALL') { 412 ss <- lsamp$sigs[,icol] 413 } else { 414 ikey <- which(lsamp$class.set == cls) 415 if (length(ikey) != 1) { 416 err.AFNI(paste("CLass", cls,"not found")) 417 return(NULL); 418 } 419 key <- lsamp$classkey.set[ikey] 420 if (verb > 1) note.AFNI(paste("Key ",key, "matches class", cls)) 421 422 ss <- lsamp$sigs[which(lsamp$labs==key),icol] 423 } 424 if (scl) { 425 sclp <- lsamp$feat.scale[feat][[1]] 426 ss <- (ss*sclp$s+sclp$d) 427 attr(ss,"scale.params") <- sclp 428 } 429 430 attr(ss,"feature") <- feat 431 attr(ss,"seg.class") <- cls 432 433 return(ss) 434} 435 436compute.feature.scaling <- function (lsamp, a=1, b=10) { 437 fs <- list() 438 if (0) { 439 for (i in 1:ncol(lsamp$sigs)) { 440 d <- min(lsamp$sigs[,i]) 441 s <- 10/(max(lsamp$sigs[,i])-d) 442 kk <- list(d=d, s= s) 443 fs <- c(fs, list(feat.scale = kk)) 444 } 445 } else { 446 for (i in 1:ncol(lsamp$sigs)) { 447 MI <- min(lsamp$sigs[,i]) 448 MX <- max(lsamp$sigs[,i]) 449 d <- (a*MX-b*MI)/(MX-MI) 450 s <- (b-a)/(MX-MI) 451 kk <- list(d=d, s= s) 452 fs <- c(fs, list(feat.scale = kk)) 453 } 454 } 455 names(fs) <- lsamp$feat.labels 456 return(fs) 457} 458 459#Write distribution parameters to niml file 460NI.write.train.dists <- function (lprob, fname='') { 461 #form the feature scale string 462 for (i in 1:length(lprob$feat.scale)) { 463 kk <- paste(sprintf('%s_Scale+Shift ="', lprob$feat.labels[i]), 464 lprob$feat.scale[i][[1]]$s, ' , ' ,lprob$feat.scale[i][[1]]$d, 465 '"\n',sep='') 466 if (i==1) tt<-kk 467 else tt <- paste(tt, kk, sep='') 468 } 469 hd <- paste('<TRAIN_DISTS\n', 470 'ni_type="2*String,2*float"\n', 471 'ni_dimen="', nrow(lprob$pfci)*ncol(lprob$pfci),'"\n', 472 'ColumnLabels="Feature ; Class ; shape ; rate"\n', 473 'ClassSet="', paste(lprob$class.set,collapse=' ; '),'"\n', 474 'FeatureSet="', paste(lprob$feat.labels,collapse=' ; '),'"\n', 475 'Dist="',lprob$pfc[1][[1]]$distname,'"\n', 476 tt, 477 '>', 478 sep='') 479 480 cat(hd, file=fname, sep='\n', append=FALSE) 481 482 fn <- rownames(lprob$pfci) 483 cn <- colnames(lprob$pfci) 484 for (f in 1:nrow(lprob$pfci)) { 485 for (c in 1:ncol(lprob$pfci)) { 486 i <- lprob$pfci[f,c] 487 par <- get.train.pdist(lprob, f, c, 0) 488 ll <- paste('"',fn[f],'" "', cn[c],'"', 489 ' ', par$estimate['shape'],' ', par$estimate['rate'], 490 sep = '') 491 cat(ll, file=fname, sep='\n', append=TRUE) 492 } 493 } 494 cat('</TRAIN_DISTS>', file=fname, sep='\n', append=TRUE) 495} 496 497 498#Get the parameters of the gamma fit for some feature, and class 499get.train.pdist <- function (lsamp, feat=NULL, cls=NULL, verb=1) { 500 pdist <- lsamp$pfc[lsamp$pfci[feat,cls]][[1]] 501 if (verb) plot(pdist) 502 return(pdist) 503} 504 505pdf.gam <- function (x,ash,brt, uns=0) { 506 if (uns) { 507 an <- (brt^ash)/gamma(ash)*(x^(ash-1))*exp(-brt*x) 508 } else { 509 lan <- ash*log(brt)-lgamma(ash)+(ash-1)*log(x)-brt*x; 510 an <- exp(lan) 511 } 512 return(an) 513} 514 515area.gam <- function (x, binw=0.1, shape, rate, meth=1) { 516 #reference 517 if (meth != 2) 518 a0 <- (pgamma(x-binw/2.0, 519 shape=shape, rate=rate, 520 lower.tail=FALSE) - 521 pgamma(x+binw/2.0, 522 shape=shape, rate=rate, 523 lower.tail=FALSE)); 524 #approx 525 if (meth != 1) 526 a1 <- (pdf.gam((x-binw/2.0), shape, rate)+pdf.gam((x+binw/2.0), shape, rate))/2.0*binw; 527 528 if (meth == 0) { 529 #debug 530 cat ('Reference, Approx = ', a0, ',', a1, '\n'); 531 return(-1.0); 532 } else if (meth == 1) { 533 #reference 534 return(a0) 535 } else { 536 #approx 537 return(a1) 538 } 539} 540 541disp.train.pdist <- function (lsamp, feats=NULL, clss=NULL, stp=0.1) { 542 if (is.null(feats)) feats <- lsamp$feat.labels 543 if (is.null(clss)) clss <- lsamp$class.set 544 imax <- length(clss)*length(feats) 545 k<-1 546 for (cls in clss) { 547 for (feat in feats) { 548 k <- k + 1 549 plot(par <- get.train.pdist(lsamp, feat, cls, verb=0)) 550 s<-0.0; 551 for (x in seq(1,10,stp)) 552 s<-s+pdf.gam(x,20.8156340614560, 3.5370637717458)*stp; 553 if (k<imax) ans <- prompt.AFNI(paste("Showing ", feat, cls, 554 "[shape=",par$estimate['shape'], 555 ',rate=' , par$estimate['rate'], 556 ', sum = ', s, 557 "]\n Hit enter to proceed.")) 558 if (ans != 1) return(0) 559 } 560 } 561} 562 563#Estimate the probability of feature amplitude given a class 564# 4-5 secs per call 565p.a_GIV_cvfu <- function (af , feat, cls, lsamp, binwidth=0.01, uid, verb=1) { 566 #note.AFNI(paste("p.a_GIV_cvfu in")); 567 fsave <- sprintf('/tmp/%s.p.a_GIV_cvfu-%s-%s.Rdat', 568 uid, feat, cls) 569 if (file.exists(fsave)) { 570 load(fsave) 571 } else { 572 par <- get.train.pdist(lsamp, feat, cls, verb=0) 573 #apply feature scaling 574 sclp <- lsamp$feat.scale[feat][[1]] 575 af <- (af*sclp$s+sclp$d) 576 p <- pgamma(af, 577 shape=par$estimate['shape'], rate=par$estimate['rate'], 578 lower.tail=FALSE) - 579 pgamma(af+binwidth, 580 shape=par$estimate['shape'], rate=par$estimate['rate'], 581 lower.tail=FALSE) 582 save(p, file=fsave) 583 } 584 #note.AFNI(paste("p.a_GIV_cvfu out")); 585 return(p) 586} 587 588#Estimate the probability of a class, given a feature 589p.cv_GIV_afu <- function (af, feat, cls, lsamp, mixfrac, uid, verb=0) { 590 ss <- 0 591 nn <- 0 592 for (v in lsamp$class.set) { 593 #note.AFNI(paste("p.cv_GIV_afu", v)) 594 pp <- mixfrac[v] * p.a_GIV_cvfu(af,feat=feat,cls=v, 595 lsamp=lsamp, uid=uid, verb=0) 596 if (v == cls) nn <- pp 597 ss <- ss + pp 598 } 599 return(nn/ss) 600} 601 602#Estimate the probability of a class, given all features 603p.cv_GIV_A <- function (A, Afeat, cls, lsamp, mixfrac, uid, verb=0) { 604 p <- 0 605 for (i in 1:ncol(A)) { 606 note.AFNI(paste("p.cv_GIV_A", i)) 607 p <- p + log(p.cv_GIV_afu(A[,i], feat=Afeat[i], 608 cls=cls, lsamp=lsamp, mixfrac=mixfrac, uid=uid, 609 verb=0)) 610 } 611 return(exp(p)) 612} 613 614disp.train.samples <- function (smp, fitprms=NULL) { 615 h<-hist(smp,breaks=100) 616 xhist<-c(min(h$breaks),h$breaks) 617 yhist<-c(0,h$density,0) 618 if (is.null(fitprms)) yfit <- NULL 619 if (is.null(atr <- attr(smp,"scale.params"))) xlab <- "feature value" 620 else { 621 xlab <- paste("feature value scaled d,s=[", atr$d, ',', atr$s, "]") 622 } 623 plot( xhist,yhist, 624 type="s",ylim=c(0,max(yhist,yfit)), 625 xlab=xlab,ylab="Density",cex.lab=1.5,cex.axis=1.5) 626 title(paste('pdf(',attr(smp,"feature"),',',attr(smp,"seg.class"),')')) 627 628} 629 630lp.norm <- function (dv) { 631 if (0) { 632 #method 0, unstable 633 dvn <- vector(length=length(dv)) 634 ddf <- 0.0 635 for (i in 1:length(dv)) { 636 ddf <- ddf+exp(dv[i]) 637 } 638 for (i in 1:length(dv)) { 639 dvn[i] = exp(dv[i])/ddf 640 } 641 } else { 642 #method 1 643 dvn <- vector(length=length(dv)) 644 for (i in 1:length(dv)) { 645 ddf <- 1.0 646 for (ii in 1:length(dv)) { 647 if (ii != i) ddf <- ddf + exp(dv[ii]-dv[i]) 648 } 649 dvn[i] = 1.0 / ddf 650 } 651 } 652 return(dvn) 653} 654 655fit.pdist <- function (smp, dist="gamma", verb=1) { 656 657 parms <- NULL 658 if (dist == "gamma") { 659 if (!exists('fitdist')) { #If we do not have fitdistrplus 660 parms<-fitdistr(as.vector(smp),"gamma", lower=0.001) 661 } else { 662 parms<-try(fitdist(as.vector(smp),"gamma")) 663 } 664 } else if (dist == "norm") { 665 if (!exists('fitdist')) { #If we do not have fitdistrplus 666 parms<-fitdistr(as.vector(smp),"norm", lower=0.001) 667 } else { 668 parms<-fitdist(as.vector(smp),"norm") 669 } 670 } else { 671 err.AFNI("Bad distribution name"); 672 return(NULL) 673 } 674 675 if (attr(parms,"class") != 'fitdist' && attr(parms,"class") != 'fitdistr' ) { 676 warn.AFNI(paste("Failed to fit sample ", 677 attr(smp,"feature"),',',attr(smp,"seg.class"))) 678 parms <- NULL 679 } 680 681 return(parms) 682 683} 684