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