.packageName <- "ppc"
ppc.cv <- function(ppc.fit, data,   user.parms){
  ## K-fold cross-validation for ppc method
  ## takes results of make.centroids.list, find.splits, and predict.ppc
  ## and computes cross-validated predictions and error estimate.
  
  peaklist.tr<- data$peaklist
  ytr <- data$ytr
  logmz <- data$logmz
  
  peak.gap <- user.parms$peak.gap
  nsplits <- user.parms$nsplits
  fix.at.one <- user.parms$fix.at.one
  recluster <- user.parms$recluster
  
  folds <- balanced.folds(ytr)
  
  n.class<-table(ytr)
  n.threshold<- length(ppc.fit$threshold)
  
  n <- length(peaklist.tr)
  yhatcv <- array(NA,c(n,n.threshold))
  probcv <- array(1, c(n, length(n.class), n.threshold))
  
  for(ii in 1:length(folds)){
    gg <- folds[[ii]]
    if (ppc.options$debug) cat(c("fold=",ii),fill=T)
    
    pk <- peaklist.tr[-gg]
    
    data.temp <- list(ytr=ytr[-gg], logmz=logmz, peaklist=pk)
    astar <- ppc.make.centroid.list(data.temp,  user.parms)
    
    junk0 <- ppc.predict.peaks(astar,  data.temp)
    aa <- ppc.find.splits(astar, junk0 ,data.temp,user.parms)
    ress <- ppc.predict(astar,aa,logmz,peaklist.te=peaklist.tr[gg],threshold=ppc.fit$threshold, metric="euclidean")
    yhatcv[gg,] <- ress$yhat
    probcv[gg,,] <- ress$prob
  }
  
  junk <- ppc.cv.error(ytr, yhatcv, folds)
  
  return(list(err=junk$err,
              se=junk$se,
              confusion=junk$confusion,
              threshold=ppc.fit$threshold,
              yhat=yhatcv,
              prob=probcv,
              y=ytr,
              folds=folds,
              numsites=ppc.fit$numsites))
}
##
## ppc.cv.error
## Compute error rate, se of error rate and confusion matrices
## from a ppc.cv fit for all values of threshold
##

ppc.cv.error<- function(y, yhat, folds) {
  temp <- ppc.error(y, yhat)
  err <- temp$err
  confusion <- temp$confusion
 
  nfolds <- length(folds)
  n <- length(y)

  codesy <- names(table(y))
  yhat2 <- matrix(codesy[yhat], nrow = nrow(yhat), ncol = ncol(yhat))
  err2 <- matrix(NA, ncol = ncol(yhat), nrow = nfolds)
  temp <- matrix(y, ncol = ncol(yhat), nrow = n)
  ni <- rep(NA, nfolds)

  for(i in 1:nfolds) {
    ii <- folds[[i]]
    ni[i] <- length(folds[[i]])
    err2[i,  ] <- apply(temp[ii,  ] != yhat2[ii,  ], 2, sum) / ni[i]
  }
  se <- sqrt(apply(err2, 2, var) / nfolds)

  return(list(err = err / n,
              confusion = confusion,
              se = se))

}


ppc.error<- function(y,yhat){
  ## compute errors and confusion matrices, from a ppc fit
  
  tt <- vector("list",ncol(yhat))
  codesy<-levels(y)
  yhat2<-matrix(codesy[yhat], nrow=nrow(yhat),ncol=ncol(yhat))
  
  
  for(i in 1:length(tt)){
    tt[[i]] <- table(y,yhat2[,i])
    
  }
  err <- rep(NA,length(tt))
  
  for(i in 1:length(tt)){
    err[i] <- sum(yhat2[,i]!=y)
  }
  
  return(list(err=err,confusion=tt))
}


##
## Change by naras. Added split.fit argument.
##
ppc.fdr <- function(data, centroid.fit, peak.fit, split.fit, ppc.fit, user.parms){
  ytr <- data$ytr
  fix.at.one <- user.parms$fix.at.one
  nsplits <- user.parms$nsplits
  nperms <- user.parms$nperms
  m <- nrow(centroid.fit$cen)
  
  threshold <- ppc.fit$threshold
  tt <- split.fit$prhat[,1]-split.fit$prhat[,2]
  
  ress <- vector("list",nperms)
  ttstar <- matrix(NA,nrow=length(tt),ncol=nperms)
  
  for(i in 1:nperms){
    
    if (ppc.options$debug) cat(c("i=",i),fill=T)
    ytr2 <- sample(ytr)
    data2 <- data
    data2$ytr <- ytr2
    split.fit2 <- ppc.find.splits(centroid.fit, peak.fit, data2, user.parms) 
    
    ress[[i]] <- split.fit2$prhat
    ttstar[,i] <- split.fit2$prhat[,1]-split.fit2$prhat[,2]
  }
  
  nt <- length(threshold)
  npeaks <- rep(NA,nt)
  npeaks0 <- npeaks
  for(i in 1:nt){
    npeaks[i] <- sum(abs(split.fit$prhat[,1]-split.fit$prhat[,2])>threshold[i])
    npeaks0[i] <- 0
    for(j in 1:nperms){
      npeaks0[i] <- npeaks0[i]+sum(abs(ress[[j]][,1]-ress[[j]][,2])>threshold[i])
    }}
  
  fdr <- (npeaks0/nperms)/npeaks
  q1 <- quantile(abs(tt), .25)
  q2 <- quantile(abs(tt), .75)
  
  pi0 <- min((sum(abs(ttstar)> q1 & abs(ttstar)< q2)/nperms)/(.5*m) ,1 )
  
  fdr <- fdr*pi0
  
  results <- cbind(threshold, npeaks,fdr)
  
  dimnames(results) <- list(NULL,c("threshold", "npeaks","fdr"))
  return(list(results=results,pi0=pi0, threshold=threshold))
}
 ppc.find.peaks  <- function(mz,x,user.parms){
##
## peak finding algorithm for SELDI spectra, from 
##  description of Yasui  et al Biostatistics 2003
##
## x is a single spectrum
##  span is window size for supersmoother estimate of background
##         minht is min height for a peak; stn is  min signal to noise
##         ratio for a peak ( i.e. must be > stn* background)
  
  span<-user.parms$span
  minht<-user.parms$minht
  
  stn<-user.parms$stn
  smoothing.span<-user.parms$smoothing.span
  
  a<-ppc.peaks(x,span=span)
  b<-supsmu(mz,x,span=smoothing.span)
  peaks2<-(1:length(mz))[a & (x> stn*b$y) & (x>minht)]
  return (cbind(mz[peaks2],x[peaks2]))
}
ppc.find.splits<- function(centroid.fit, peak.fit, data, user.parms)
{
  ## find best split points for training data
  ## takes centroids.fit- result of call to make.centroids.list
  ## and peaks.fit- result of call to predict.peaks
  ##
  ##  ytr is  the vector of  class labels 
  ## nsplits is number of equally spaced  split points to try (plus the value 0)
  ##  fix.at.one=-TRUE means zero is the only split point tried
  ##    (ie no peak vs peak)
  ##
  ## returns prhat - estimated optimal split proportions for each site, 
  ##  pr- proportions in each  split category for each class, 
  ## prclose- logical matrix indicating split values that coem within 10% of best
  ## cuthat,cutpoints- optimal cut ## and matrix of all cutpoints considered
  ##
  
  ytr <- data$ytr
  
  n.class<- table(ytr)
  
  nsplits <- user.parms$nsplits
  fix.at.one <- user.parms$fix.at.one 
  
  which.is.max.na <- function(x)
    {
      xx <- x[!is.na(x)]
      y <- seq(length(xx))[xx == max(xx)]
      if(length(y) > 1) {
        y <- sample(y, 1)
      }
      o <- (1:length(x))[!is.na(x)]
      return(o[y])
    }
  ht <- peak.fit$ht * peak.fit$ind
  cent <- centroid.fit$cent
  alpha <- (1:nsplits)/(nsplits + 1)
  p <- nrow(ht)
  Y <- model.matrix( ~ factor(ytr) - 1, data = list(ytr = ytr))
  K <- length(table(ytr))
  pr <- array(0, c(p, nsplits + 1, K))
  qu <- matrix(NA, nrow = p, ncol = nsplits + 1)
  
  ##at each site, compute proportion  of peak heights exceeding the split  quantiles 
  
  for(i in 1:p) {
    if (ppc.options$debug) cat(i,fill=T)
    o <- ht[i,  ] > 0
    temp <- sort(ht[i, o])
    temp2 <- pmax(1, trunc(length(temp) * alpha))
    qu[i,  ] <- c(0, temp[temp2])
    for(j in 1:(nsplits + 1)) {
      for(k in 1:K) {
        pr[i, j, k] <- sum(ht[i, Y[, k] == 1] > qu[i, j])/sum(Y[, k] == 1)
      }
    }
  }
  prmean <- apply(pr, c(1, 2), mean)
  
  ## find best split indices
  
  if(!fix.at.one) {
    prtemp <- apply(abs(pr - array(prmean, c(p, nsplits + 1, K))),
                    c(1, 2), sum)
    cuthat <- apply(prtemp, 1, which.is.max.na)
  }
  else {
    cuthat <- rep(1, p)
  }
  prhat <- matrix(NA, nrow = p, ncol = K)
  
  
  ##compute optimized proportions and  splits within 10% of the best
  
  if(fix.at.one){prclose <- NULL}
  
  
  prclose <- matrix(F, nrow = p, ncol = nsplits + 1)
  for(i in 1:p) {
    prhat[i,  ] <- pr[i, cuthat[i],  ]
    if(!fix.at.one){
      for(j in 1:(nsplits + 1)) {
        prclose[i, j] <- prtemp[i, j] > 0.9 * prtemp[i, cuthat[
                                                               i]]
      }}
  }
  
  return(list(prhat = prhat, pr = pr, n.class=n.class, cutpoints = qu, cuthat = cuthat, prclose
              = prclose,nsplits=nsplits, fix.at.one=fix.at.one))
}


ppc.make.centroid.list<-
function (data, user.parms) 
{
    peaklist <- data$peaklist
    ytr <- data$ytr
    logmz <- data$logmz
    recluster <- user.parms$recluster
    peak.gap <- user.parms$peak.gap
    n <- length(ytr)
    m <- NULL
    for (i in 1:n) {
        m <- c(m, peaklist[[i]][, 1])
    }
    mm <- sort(unique(m))
    clust.tree <- hclust.1d(mm, debug = as.integer(ppc.options$debug))
    if (ppc.options$debug) 
        cat("bef ppc.make.centroids", fill = TRUE)
    cent <- ppc.make.centroids(clust.tree, mm, user.parms)
    if (ppc.options$debug) 
        cat("aft ppc.make.centroids", fill = TRUE)
    return(list(cent = cent, peaklist = peaklist, clust.tree = clust.tree, 
        all.peaks = mm, peak.gap = user.parms$peak.gap, recluster = user.parms$recluster))
}
ppc.make.centroids <-function(clust.tree, x, user.parms)
{
  peak.gap <-  user.parms$peak.gap
  nclust <- user.parms$nclust
  recluster <- user.parms$recluster
  
  if(!is.null(peak.gap) & !is.null(nclust))
    stop("Error: only one of peak.gap  and nlcust should be specified"
         )
  if(!is.null(nclust)) {
    aa <- cutree(clust.tree, k = nclust)
  }
  if(!is.null(peak.gap)) {
    aa <- cutree(clust.tree, h = peak.gap)
    n <- nrow(clust.tree$merge) + 1
    uaa <- unique(aa)
    cen <- matrix(NA, ncol = 4, nrow = length(uaa))
    i <- 0
    for(ii in uaa) {
      b <- x[(1:n)[aa == ii]]
      i <- i + 1
                                        #     cen[i, 1] <- (min(b) + max(b))/2
                                        #      cen[i, 2] <- (max(b) - min(b))/2
      if (ppc.options$debug) cat(i, fill = T)
      cen[i, 1] <- medoid(b)
                                        # i changed this on oct 27
                                        # cen[i, 2] <- max(cen[i,1] - min(b), max(b) - cen[i,1])
      cen[i, 2] <- max(b) - min(b)
      cen[i, 3] <- min(b)
      cen[i,4] <-  max(b)
    }
  }
  
  if(recluster & !is.null(peak.gap)){
    mind<- min(diff(cen[,1]))
    while(mind < peak.gap){
      if (ppc.options$debug) cat("reclustering",fill=T)
      xx <- cen
      clust.tree<- hclust.1d(xx[,1], debug=ppc.options$debug)
      aa <- cutree(clust.tree, h = peak.gap)
      n <- nrow(clust.tree$merge) + 1
      uaa <- unique(aa)
      cen <- matrix(NA, ncol = 4, nrow = length(uaa))
      i <- 0
      for(ii in uaa) {
        b <- xx[,1][(1:n)[aa == ii]]
        bmin <- xx[,3][(1:n)[aa == ii]]
        bmax <- xx[,4][(1:n)[aa == ii]]
        i <- i+1
        
        cen[i, 1] <- mean(b)
        cen[i, 3] <- min(bmin)
        cen[i,4] <-  max(bmax)
        cen[i,2] <- cen[i,4]-cen[i,3]
      }
      mind<- min(diff(cen[,1]))
    }
  }
  
  dimnames(cen) <- list(NULL, c("pos", "diameter", "min", "max"))
  o <- order(cen[, 1])
  return(cen[o,  ])
}

ppc.make.peaklist <- function(data, user.parms){
  xtr <- data$xtr
  mz <- data$mz
  logmz<- data$logmz
  
  pk1 <- vector("list",ncol(xtr))
  ii <- 0
  
  for(i in 1:ncol(xtr)) {
    a <- ppc.find.peaks(mz,xtr[,i],user.parms)
    aa<-match(a[,1],mz)
    ii <- ii+1
    if (ppc.options$debug) cat(ii)
    pk1[[ii]] <- cbind(logmz[aa],a[,2])
  }
  
  return(pk1)
}
##
## This function min.na is not needed as the function min offers a direct
## facility for not considering NAs.
##
##min.na<-function(x) {
##  min(x[!is.na(x)])
##}

which.is.min <- function(x) {
  y <- seq(length(x))[x == min(x, na.rm=TRUE)]
  y <- y[!is.na(y)]
  if(length(y) > 1)
    y <- sample(y, 1)
  y
}

medoid <- function(x) {
  n <- length(x)
  if(n==1){med <- x}
  if(n>1){
    d <- dist(x)
    dd <- matrix(0,nrow=n,ncol=n)
    dd[row(dd)>col(dd)] <- d
    dd <- dd+t(dd)
    
    m <- apply(dd,2,sum)
    med <- x[(1:n)[m==min(m)]]
    if(length(med)>1){
      med <- mean(med)
    }
  }
  return(med)
}

permute.rows <-function(x) {
  dd <- dim(x)
  n <- dd[1]
  p <- dd[2]
  mm <- runif(length(x)) + rep(seq(n) * 10, rep(p, n))
  matrix(t(x)[order(mm)], n, p, byrow = TRUE)
}


balanced.folds <- function(y, nfolds = min(min(table(y)), 10)) {
  totals <- table(y)
  fmax <- max(totals)
  nfolds <- min(nfolds, fmax)     
  ## makes no sense to have more folds than the max class size
  folds <- as.list(seq(nfolds))
  yids <- split(seq(y), y)        
  ## nice we to get the ids in a list, split by class
  ##Make a big matrix, with enough rows to get in all the folds per class
  bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
  for(i in seq(totals)) {
    bigmat[seq(totals[i]), i] <- sample(yids[[i]])
  }
  smallmat <- matrix(bigmat, nrow = nfolds)       # reshape the matrix
  ## Now do a clever sort to mix up the NAs
  smallmat <- permute.rows(t(smallmat))
  ## Now a clever unlisting
  ## the "clever" unlist doesn't work when there are no NAs
  ##       apply(smallmat, 2, function(x)
  ##        x[!is.na(x)])
  res <-vector("list", nfolds)
  for(j in 1:nfolds) {
    jj <- !is.na(smallmat[, j])
    res[[j]] <- smallmat[jj, j]
  }
  return(res)
}

ppc.peaks <- function(x,span){
  
  
# note-  changed this so that span is now a percentage rather than a
#   number of points
 
   ispan<-trunc(length(x)*span)

  if(ispan%%2==0){ispan <- ispan+1}
  
  n <- length(x) 
  
#  dyn.load("/home/tibs/PAPERS/protein/ppc/work/PEAKS.so")
  junk <- .Fortran("peaks",
                   x,
                   as.integer(ispan),
                   as.integer(n),
                   ans=integer(n),
                   PACKAGE="ppc")
  return(junk$ans==1)
}

hclust.1d <- function(x, debug=FALSE) {
  ##(fast) complete linkage hierarhical clustering, in one dimension
  ## R. Tibshirani Nov 2003
  # to make it work under Linux,
  # rob uncommented dyn.load line, commented two lines and added scrat2=
  
  n<-length(x)
  storage.mode(x)<- "double"
  storage.mode(n)<- "integer"
  ##
  ## dyn.load not needed when packaged!
    dyn.load("/home/tibs/PAPERS/protein/ppc/work/hclust1d.so")
  ##

  idebug <- ifelse(debug, 1, 0)
  
  junk<- .Fortran("hclust1d",
                  as.integer(idebug),
                  x,
                  n,
                  merge=integer( (n-1)*2),
                  height=double(n-1),
                  order=integer(n),
                  scrat=double(n*3),
                  scrat2=double(n*3),
                  PACKAGE="ppc")

  
  merge<-matrix(junk$merge,ncol=2,byrow=F)
  order<-junk$order
  height<-junk$height
  return(list(merge=merge,height=height,order=order))
}





ppc.peak.summary <- function(centroid.fit, peak.fit, data, user.parms, split.fit) {
  temp<-  peak.fit$ind*peak.fit$ht
  temp[temp==0] <- NA
  
  sum.na<- function(x){sum(x[!is.na(x)])}
  
  splits <- rep(NA,length(split.fit$cuthat))
  for(i in 1:length(split.fit$cuthat)){
    splits[i]<- split.fit$cutpoints[i,split.fit$cuthat[i]]
  }
 nc<-ncol(split.fit$prhat)

  labs0 <- rep(NA,nc)
  for(i in 1:nc){
    labs0[i] <- paste("Prop.in.class.",as.character(i),sep="")
}
  prdiff<-matrix(NA,ncol=nc*(nc-1)/2,nrow=nrow(split.fit$prhat))
  ii<-0
  labs <- rep(NA,nc*(nc-1)/2)
  for(i in 1:(nc-1)){
    for(j in (i+1):nc){
      ii <- ii+1
  prdiff[,ii] <- split.fit$prhat[,j]-split.fit$prhat[,i]
      labs[ii] <- paste("Pr",as.character(j),"-","Pr",as.character(i),sep="")
    }}
      
  rank.split <- rank(-apply(abs(prdiff),1,sum))
  
  res <- cbind(exp(centroid.fit$cent[,-2]), apply(temp>0,1,sum.na), rank.split,splits,split.fit$prhat, prdiff, temp)
  
  dimnames(res) <- list(NULL,c("peak.position",
                               "min",
                               "max",
                               "number.of.spectra",
                               "rank.of.peak",
                               "height.split.point",
                               labs0,
                               labs,
                               data$sample.labels))
  return(res)
}



ppc.plot.hist <-  function(peak.fit, ppc.fit, centroid.fit, split.fit, data, first.site=1, last.site=25, title.plot=NULL){
# first.site and last.site  indicate how many peak  sites to plot.   default is 
#  1 thru 25

require(class)
sitelist <- ppc.fit$sites[[1]]

if(!is.null(first.site) & is.null(last.site)){
  last.site<- first.site+24
}

if(is.null(first.site) & !is.null(last.site))
             {stop(" Error: first.site missing")}


if(last.site<first.site){stop(" Error: last.site<first.site")}

if(last.site-first.site>24){ last.site<- first.site+24}

if(!is.null(first.site)){
  sitelist<-sitelist[first.site: last.site]
}

ytr<- data$ytr
codesy <- levels(ytr)

nc<- dim(split.fit$pr)[3]
nrows <- trunc(20/nc)
par(mfcol=c(nrows,5))

if(!is.null(title.plot)){ par(oma=c(0,0,5,0))}
par(mar=c(0,1,0,1))
par(cex=.6)
ht <- peak.fit$ind*peak.fit$ht


for(i in sitelist){
 
  if (ppc.options$debug) cat(i,fill=T)
  xmax <- max(ht[i,])

  ## Ensuring that the breaks in each histogram are same for all classes
  br <- hist(ht[i,],col=3,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=F, main="")$br
  
  par(mar=c(0,1,1,1))
  junk <- hist(ht[i,ytr==codesy[1]],col=3,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=F
               ,breaks=br)
  hist(ht[i,ytr==codesy[1]],col=5,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br, main="")
  box(col=gray)
  wid <- (junk$br[2]-junk$br[1])
  v <- split.fit$cutpoints[i,split.fit$cuthat[i]]+wid
  v0<-v
  if(v==0){v<-v+wid}
  
  h <- max(junk$counts)
  
  abline(v=v,lty=1,col=2)
  
  for(j in (1:ncol(split.fit$prclose))[split.fit$prclose[i,]]){
    points(split.fit$cutpoints[i,j]+wid,h/2,col=2,pch="x")
  }
  
  ## Pattern match to figure 
  o<-knn1(data$logmz, centroid.fit$cent[i,1], 1:length(data$logmz))
  lab<- data$mz.labels[o]
  text(xmax*.5,.9*h,labels=lab,cex=.8)
  text(v,3*h/4,labels=as.character(round(v0,2)),col=6)
  
  
  pr <- sum(ht[i,ytr==codesy[1]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[1])
  text(xmax,h/2,labels=as.character(round(pr,2)),col=1,cex=.7)
  
  if(nc>2){
    for(ii in 2:(nc-1)){
      par(mar=c(1,1,1,1))
      junk2 <- hist(ht[i,ytr==codesy[ii]],col=4,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=
                    F,breaks=br)
      hist(ht[i,ytr==codesy[ii]],col=ii,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br, main="") 
      
      box(col=gray)
      abline(v=v,lty=1,col=2)
      h2 <- max(junk2$counts)
      pr <- sum(ht[i,ytr==codesy[ii]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[ii])
      text(xmax,h2/2,labels=as.character(round(pr,2)),col=1,cex=.7)
}}

par(mar=c(2,1,0,1))
junk2 <- hist(ht[i,ytr==codesy[nc]],col=4,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=
F,breaks=br)
hist(ht[i,ytr==codesy[nc]],col=nc,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br,main="") 
  
box(col=gray)
abline(v=v,lty=1,col=2)
h2 <- max(junk2$counts)
  pr <- sum(ht[i,ytr==codesy[nc]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[nc])
text(xmax,h2/2,labels=as.character(round(pr,2)),col=1,cex=.7)

}

if(!is.null(title.plot)){ mtext(title.plot,outer=TRUE,side=3)}

return()
}

ppc.plotcv <- function(fit) {
  par(mar = c(5, 5, 5, 1))
  par(mfrow = c(2, 1))
  n <- nrow(fit$yhat)
  y <- fit$y
  nc <- length(table(y))
  nfolds <- length(fit$folds)

  plot(fit$threshold, fit$err, ylim = c(-0.1, 0.8), xlab = 
       "Value of threshold  ", ylab = "Misclassification Error", type
       = "n", yaxt = "n")
  axis(3, at = fit$threshold, lab = paste(fit$numsites), srt = 90, adj = 0)
  mtext("Number of sites", 3, 4, cex = 1.2)
  axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
  lines(fit$threshold, fit$err, col = 2)
  o <- fit$err == min(fit$err)
  points(fit$threshold[o], fit$err[o], pch = "x")
  error.bars(fit$threshold, fit$err - fit$se, fit$err + fit$se)
  err2 <- matrix(NA, nrow = length(unique(y)), ncol = length(fit$threshold
                                                 ))
  for(i in 1:(length(fit$threshold) - 1)) {
    s <- fit$confusion[[i]]
    diag(s) <- 0
    err2[, i] <- apply(s, 1, sum)/table(y)
  }
  plot(fit$threshold, err2[1,  ], ylim = c(-0.1, 1.1), xlab = 
       "Value of threshold ", ylab = "Misclassification Error", type
       = "n", yaxt = "n")
  axis(3, at = fit$threshold, lab = paste(fit$numsites), srt = 90, adj = 0)   
 # mtext("Number of sites", 3, 4,cex=1.2)
  axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
  for(i in 1:nrow(err2)) {
    lines(fit$threshold, err2[i,  ], col = i + 1)
  }
  legend(0, 0.9, dimnames(table(y))[[1]], col = (2:(nc + 1)), lty = 1)
  par(mfrow = c(1, 1))
}

error.bars <-function(x, upper, lower, width = 0.02, ...) {
  xlim <- range(x)
  barw <- diff(xlim) * width
  segments(x, upper, x, lower, ...)
  segments(x - barw, upper, x + barw, upper, ...)
  segments(x - barw, lower, x + barw, lower, ...)
  range(upper, lower)
}

ppc.plotcvprob <- function(fit, data, threshold) {
  par(pch = 1)
  ii <- (1:length(fit$threshold))[fit$threshold > threshold]
  ii <- ii[1]
  ss <- data$samplelabels
  pp <- fit$prob[,  , ii]

  y <- fit$y

  o <- order(y)
  y <- y[o]
  if(!is.null(ss)) {
    ss <- ss[o]
  }
  ppp <- pp[o,  ]
  n <- nrow(ppp)
  nc <- length(unique(y))
  par(cex = 1)
  plot(1:n, ppp[, 2], type = "n", xlab = "sample", ylab = 
       "cross-validated probabilities",  axes = FALSE)
  axis(1)
  labs <- round(seq(min(ppp),max(ppp),length=5),2)
  axis(2, labels = as.character(labs),at=labs)
  for(j in 1:nc) {
    points(1:n, ppp[, j], col = j + 1)
  }
  for(j in 1:(nc - 1)) {
    abline(v = cumsum(table(y))[j] + 0.5, lty = 2)
  }
  h <- c(0, table(y))
  for(j in 2:(nc + 1)) {
    text(sum(h[1:(j - 1)]) + 0.5 * h[j], 1.02, label = levels(y)[j - 
                                                 1], col = j)
  }
  abline(h = 1)
  if(!is.null(ss)) {
    text(1:length(ss), 1.1, labels = ss, srt = 90, cex = 0.7)
  }
  ##if(!is.null(ss)){axis(3,labels=ss,at=1:length(ss),srt=90)}
}


ppc.plotfdr <- function(fdrfit){
  plot(fdrfit$results[,"npeaks"],fdrfit$results[,"fdr"],
xlab="Number of peaks called significant",
ylab="False discovery rate",type="b",log="xy")
  axis(3,at=fdrfit$results[,"npeaks"], labels=round(fdrfit$threshold,2))
  
  return()
}
ppc.predict <- function(centroid.fit, split.fit, logmz,  peaklist.te, n.threshold = 30, threshold=NULL, metric = 
                        c("binomial","euclidean", "absolute"), summ=c("mean","median"))
{
  ## test set prediction for PPC method
  ## makes  predictions  for test set list of peaks in peaklist.te
  ##  
  ## "metric" is the metric used in the nearest centroid rule; "summ" is
  ## summary used  to combine distances over sites
  ## n.threshold is the number of shrinkage thresholds used
  
  ## returns yhat, plus ind0= indicator matrix  of whether peak was found,
  ##           ind=indicator matrix  of whether peak > cutpoint was found
  ##           ht=matrix of peak heights
  
  ## "numsites" is the number of sites present at each threshold; "sites" is a list
  ## containing the ##s of these sites
  ##
  ## dis is the distance of each test profile to each class centroid
  
  ##

SMALL<- 10e-6

  metric <- match.arg(metric)
  print(metric)
summ <- match.arg(summ)
  print(summ)
  if(!is.null(threshold)){n.thresholds<-length(threshold)}
  
  n <- length(peaklist.te)
  n.class <- split.fit$n.class
  K<- ncol(split.fit$prhat)
  p<-length(split.fit$cuthat)
  yhat <- matrix(NA, nrow = p, ncol = n)
  yhatt <- matrix(NA, ncol = n.threshold, nrow = n)
  ind0 <- yhat
  ht <- yhat
  dis <- array(NA,c(n, K, n.threshold))
  
  numsites <- rep(NA, n.threshold)
  sites <- vector("list", n.threshold)
  for(j in 1:n) {
    if (ppc.options$debug) cat(j)
    aaa <- ppc.predict.peaks1(centroid.fit, logmz,  peaklist= 
                          peaklist.te[[j]])
    ht.hat <-  aaa$ht * aaa$ind
    ind0[,j] <- 1*aaa$ind
    ht[,j] <- ht.hat
    for(i in 1:p) {
      yhat[i, j] <- 1 * (ht.hat[i] > split.fit$cutpoints[i, split.fit$cuthat[i]])
    }
    
  }        
  
  prmean <- apply(split.fit$prhat,1,mean)
  delta <- split.fit$prhat-matrix(prmean,nrow=nrow(split.fit$prhat),ncol=K)
  
  dd<- matrix(NA,nrow=n,ncol=K)
  
  if(is.null(threshold)){threshold <- seq(0, max(abs(delta)), length = n.threshold)}
  
  soft.thresh <- function(x, tt)
    {
      sign(x) * (abs(x) - tt) * (abs(x) > tt)
    }
  
  for(j in 1:n.threshold) {
    if (ppc.options$debug) cat(j)
    delta2 <- soft.thresh(delta, threshold[j])
    pr <- prmean + delta2
    sumabs<-  apply(abs(delta2),1,sum)
    pos <-  sumabs> SMALL
    numsites[j] <- sum(pos)
    sites[[j]] <- (1:nrow(pr))[pos]
    
    sites[[j]] <- sites[[j]][order( - abs(sumabs[pos]))]
    if(metric == "euclidean") {
      for(k in 1:K){
        temp<- (yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))^2
        dd[,k] <- apply(temp,2,summ)
      }
      
    }
    if(metric == "absolute") {
      for(k in 1:K){
        temp<- abs(yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))
        dd[,k] <- apply(temp,2,summ)
      }
    }
    if(metric == "binomial") {
      
      ## compute Bayes estimate of pr, to keep it away from 0 or 1
      ntemp <- matrix(n.class,nrow=nrow(pr),ncol=length(n.class), byrow=T)
      
      prbayes<- (pr*ntemp+1)/(ntemp+2)
      
      for(k in 1:K){
        prmat<-matrix(prbayes[,k], ncol = n, nrow = nrow(yhat))
        temp<- -1*(yhat*log(prmat)+(1-yhat)*log(1-prmat))
        dd[,k] <- apply(temp,2,summ)
      }
      
    }
    
    yhatt[, j] <- apply(dd,1,which.is.min)
    dis[,,j] <- dd
    
  }
  prob <- exp(-dis)
  prob <- prob/array(apply(prob,c(1,3),sum),dim(prob))
  
  return(list(yhat = yhatt, threshold = threshold, numsites = numsites, sites
              = sites, ind = yhat, ind0=ind0, prob=prob, ht=ht))
  
}
ppc.predict.peaks <- function(centroid.fit, data)
{
  ## takes centroid.fit (result of a call to make.centroids.list)
  ##  and looks for these m peaks  in peaklist, a list of length n
  ## returns   ind - an m by n matrix of TRUE/FALSE values
  ##       and ht, the matrix of  corresponding peak heights
  
  
peaklist<-data$peaklist
logmz<- data$logmz

  res2 <- NULL
  ht <- NULL
  for(i in 1:length(peaklist)) {
    if (ppc.options$debug) cat(i)
    
    junk <- ppc.predict.peaks1(centroid.fit, logmz,  peaklist.new = 
                           peaklist[[i]])
    
    res2 <- cbind(res2, junk$ind)
    ht <- cbind(ht, junk$ht)
    
  }
  return(list(ind = res2, ht = ht))
}
ppc.predict.peaks1 <-  function(centroid.fit, logmz,  peaklist.new)
{
  require(class)
  ## looks in peaklist.new for peaks that are near the centroid peaks
  ##   defined in centroid.fit
  ## returns logical indicator vector "ind" and peak heights "ht"
  
# changed to work with new test data with new mz values
#  pee <- match(peaklist.new[, 1], logmz)

pee<-knn1(matrix(logmz,ncol=1),peaklist.new[, 1,drop=FALSE], 1:length(logmz))
  
  xnew <- rep(NA, length(logmz))
  xnew[pee] <- peaklist.new[, 2]
  
  m3 <- logmz[pee]
  cent <- centroid.fit$cent
 
  bb <- knn1(matrix(m3, ncol = 1), cent[, 1, drop = F], 1:length(m3))
  
  ind <- abs(m3[bb] - cent[, 1]) <=centroid.fit$peak.gap/2 
  ht <- xnew[pee][bb]

#NOTE: ht is returned as non-zero, even if ind=0 (ie peak does not match
#a training peak)

  return(list( ind =ind, ht = ht))
}
ppc.predict1 <- function(centroid.fit, split.fit, logmz,  peaklist.te,  threshold, metric = 
                        c("binomial","euclidean", "absolute"), summ=c("mean","median"))
{
  ## test set prediction for PPC method
  ## makes  predictions  for test set list of peaks in peaklist.te
  ##  make predictions for a single value of threshold, while
  #    ppc.predict makes prediction for a set of thresholds
  ##  
  ## "metric" is the metric used in the nearest centroid rule; "summ" is
  ## summary used  to combine distances over sites
  ## n.threshold is the number of shrinkage thresholds used
  
  ## returns yhat, plus ind0= indicator matrix  of whether peak was found,
  ##           ind=indicator matrix  of whether peak > cutpoint was found
  ##           ht=matrix of peak heights
  
  ## "numsites" is the number of sites present at each threshold; "sites" is a list
  ## containing the ##s of these sites
  ##
  ## dis is the distance of each test profile to each class centroid
  
  ##
  metric <- match.arg(metric)
  summ <- match.arg(summ)
  
  if(!is.null(threshold)){n.thresholds<-length(threshold)}
  
  n <- length(peaklist.te)
  n.class <- split.fit$n.class
  K<- ncol(split.fit$prhat)
  p<-length(split.fit$cuthat)
  yhat <- matrix(NA, nrow = p, ncol = n)
  yhatt <- rep(NA,  n)
  ind0 <- yhat
  ht <- yhat
  dis <- array(NA,c(n, K))
  
  for(j in 1:n) {
#    if (ppc.options$debug) cat(j)
    aaa <- ppc.predict.peaks1(centroid.fit, logmz,  peaklist= 
                          peaklist.te[[j]])
    ht.hat <-  aaa$ht * aaa$ind
    ind0[,j] <- 1*aaa$ind
    ht[,j] <- ht.hat
    for(i in 1:p) {
      yhat[i, j] <- 1 * (ht.hat[i] > split.fit$cutpoints[i, split.fit$cuthat[i]])
    }
    
  }        
  
  prmean <- apply(split.fit$prhat,1,mean)
  delta <- split.fit$prhat-matrix(prmean,nrow=nrow(split.fit$prhat),ncol=K)
  
  dd<- matrix(NA,nrow=n,ncol=K)
  
  
  soft.thresh <- function(x, tt)
    {
      sign(x) * (abs(x) - tt) * (abs(x) > tt)
    }
  
    delta2 <- soft.thresh(delta, threshold)
    pr <- prmean + delta2
    sumabs<-  apply(abs(delta2),1,sum)
    pos <-  sumabs!=0
    numsites <- sum(pos)
    sites <- (1:nrow(pr))[pos]
    
    sites <- sites[order( - abs(sumabs[pos]))]
    if(metric == "euclidean") {
      for(k in 1:K){
        temp<- (yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))^2
        dd[,k] <- apply(temp,2,summ)
      }
      
    }
    if(metric == "absolute") {
      for(k in 1:K){
        temp<- abs(yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))
        dd[,k] <- apply(temp,2,summ)
      }
    }
    if(metric == "binomial") {
      
      ## compute Bayes estimate of pr, to keep it away from 0 or 1
      ntemp <- matrix(n.class,nrow=nrow(pr),ncol=length(n.class), byrow=T)
      
      prbayes<- (pr*ntemp+1)/(ntemp+2)
      
      for(k in 1:K){
        prmat<-matrix(prbayes[,k], ncol = n, nrow = nrow(yhat))
        temp<- -1*(yhat*log(prmat)+(1-yhat)*log(1-prmat))
        dd[,k] <- apply(temp,2,summ)
      }
      
    }
    
    yhatt <- apply(dd,1,which.is.min)
    dis <- dd
    

  prob <- exp(-dis)
  prob <- prob/array(apply(prob,1,sum),dim(prob))
  
  return(list(yhat = yhatt, threshold = threshold, numsites = numsites, sites
              = sites, ind = yhat, ind0=ind0, prob=prob, ht=ht))
  
}
ppc.read.peaks.batch<- function(dir, batches){

comm<- paste("ls ",dir,sep="")

f1<-NULL
for(i in batches){
  temp<- paste(comm, "/",i,sep="")
  temp2<-system(temp,intern=TRUE )
  temp3<-paste(dir,"/",i,"/",temp2,sep="")
  f1<-c(f1,temp3)
}


n<-length(f1)
peaklist<- vector("list",n)

logmz<-NULL
 for(i  in 1:length(f1)){
        peaklist[[i]]<- ppc.read.peaks.file(f1[i])
        logmz<-c(logmz,peaklist[[i]][,1])
}

return(list(peaklist=peaklist,filenames=f1, logmz=logmz))

}
ppc.read.peaks.file<- function(filename){
# note: we ignore anything wiht mz <=0 
 pat<-scan(filename,what="", comment.char="#")
        pat<-matrix(pat,ncol=7,byrow=T)[,(1:2)]
        pat<-matrix(as.numeric(as.character(pat)),ncol=2)
        o<-pat[,1]>0
        pat<-pat[o,]
        return( cbind(log(pat[,1]),pat[,2]))
}

ppc.read.peaks.nobatch<- function(dir){
comm<- paste("ls ",dir,sep="")
f1<-system(comm,intern=T)

n<-length(f1)
peaklist<- vector("list",n)

logmz<-NULL
 for(i  in 1:length(f1)){
        peaklist[[i]]<- ppc.read.peaks.file(paste(dir, "/",f1[i],sep=""))
  logmz<-c(logmz,peaklist[[i]][,1])
}

return(list(peaklist=peaklist,filenames=f1, logmz=logmz))

}
ppc.read.raw.batch<- function(dir, batches, mz=NULL){

comm<- paste("ls ",dir,sep="")

f1<-NULL
for(i in batches){
  temp<- paste(comm, "/",i,sep="")
  temp2<-system(temp,intern=TRUE )
  temp3<-paste(dir,"/",i,"/",temp2,sep="")
  f1<-c(f1,temp3)
}



if(is.null(mz)){
	pat1 <- read.table(f1[1],sep=",")
        mz <- pat1[,1]
}

xtr <- matrix(NA,nrow=length(mz),ncol=length(f1))
for(j in 1:length(f1)){
  if (ppc.options$debug) cat(j,fill=T)
  temp<-  read.table(f1[j],sep=",")

  xtr[,j] <-approx(temp[,1], temp[,2],xout=mz)$y
}

return(list(xtr=xtr,mz=mz,filenames=f1))
}
##
## Read raw data from a directory and return a list containing
## -- a matrix of x
## -- a vector of the mz values
## -- the list of file names
##

ppc.read.raw.nobatch <- function(directory, mz = NULL) {
  datafiles.list <- ppc.xl.get.names.of.files(directory)
  
  ##
  ## First determine the dimensions of the data matrix
  ##

  if (is.null(mz)) {
    pat1 <- read.table(datafiles.list[1])
    mz <- pat1[,1]
  }
  
  xtr <- matrix(NA, nrow=length(mz), ncol=length(datafiles.list))

  ##
  ## Now read in all the data
  ##
  for (j in 1:length(datafiles.list)) {
    if (ppc.options$debug) print(paste("Reading file", datafiles.list[j]))
    
    temp <-  read.table(datafiles.list[j], sep = ",")
    
    xtr[, j] <- approx(temp[, 1], temp[, 2], xout = mz)$y
  }
  
  return(list(xtr = xtr,
              mz = mz,
              filenames=datafiles.list))
}

ppc.remove.beforeslash.and.suffix<-
 function(x){
# this function  removes everything before the last "/" at the
#beginning of a filename,  and ".xxx" from the end of a filename
#  of a vector of filenames

for(ii in 1:length(x)){
 for(i in 1:nchar(x[ii])){
     if(substring(x[ii],i,i)=="/") val<-i
 }
 x[ii]<-substring(x[ii],val+1, nchar(x[ii]))
i<-nchar(x[ii])
 while(substring(x[ii],i,i)!="."){
    i<-i-1
     }
 x[ii]<-substring(x[ii],1,i-1)
}
return(x)
}

# Naras's rewrite of Rob's functrion. Doesn't work correctly on linux

##pcc.remove.beforeslash.and.suffix <- function(x) {
##
## this function  removes everything before the last "/" at the
## beginning of a filename,  and ".xxx" from the end of a filename
## of a vector of filenames
##
##base.names <- basename(x)
##return(sapply(base.names, ppc.remove.suffix))
##}




ppc.remove.suffix<-
 function(x){
# this function just removes a .xxx suffix from the end of a filename
#  of a vector of filenames

for(ii in 1:length(x)){
 i<-nchar(x[ii])
while(substring(x[ii],i,i)!= "."){i<-i-1}
 x[ii]<-substring(x[ii],1,i-1)
}
return(x)
}



##
## Original version
##

#ppc.subset.and.reshape<- function(data, user.parms){
  
#  xtr<-data$xtr
#  mz<-data$mz
#  logmz<-data$logmz
#  ytr<-data$ytr
  
  
#  mz.labels<-NULL
#  sample.labels <- data$sample.labels
  
#  if (!is.null(user.parms$mz.min) | !is.null(user.parms$mz.max)) {
#    o1<- user.parms$mz.min
#    o2<- user.parms$mz.max
#    if(!is.null(o1)) {
#      xtr<- xtr[mz>= o1,]
#      logmz<- logmz[mz>= o1]
#      mz<- mz[mz>= o1]
#    }
#    if (!is.null(o2)) {
#      xtr<- xtr[mz<= o2,]
#      logmz<- logmz[mz<= o2]
#      mz<-    mz[mz<= o2]
#    }
#  }
#  ## here I reshape the data into one column per patient. mz and logmz values are
#  ## string out into one long vector, spaced apart

#  if (!is.null(data$batch.labels)) {
#    npatients<-length(unique(data$patient.labels))
#    nbatches<-length(unique(data$batch.labels))
#    p<-nrow(xtr)
    
#    labelord<-order(data$batch.labels)
#    xtr<-xtr[,labelord]
#    if(!is.null(ytr)){ ytr<- ytr[labelord][1:npatients]}
    
#    xtr3<-matrix(NA,nrow=nbatches*nrow(xtr),ncol=npatients)
#    ii<-1
#    for (i in 1:nbatches) {
#      xtr3[ii:(ii+p-1),]<- xtr[, ((i-1)*npatients+1):(i*npatients)]
#      ii<-(ii+p)
#    }
    
    
#    mz4<-NULL
#    logmz<-NULL
#    mz.labels<-NULL
#    fac<- 2*(max(mz)-min(mz))
#    faclog<- 2*(max(log(mz))-min(log(mz)))
    
#    for(ii in 1:nbatches) {
#      mz4<-c(mz4,mz+(ii-1)*fac)
#      logmz<-c(logmz,log(mz)+(ii-1)*faclog)
#      mz.labels<-c(mz.labels,paste(round(mz,1),"-",batches[ii],sep=""))
#    }
#    xtr<-xtr3
#    mz<-mz4
    
#    sample.labels<-data$patient.labels[labelord]
#    sample.labels<-sample.labels[1:npatients]
#  }
  
#  iytr<-NULL
#  if(!is.null(ytr)) {
#    iytr<- as.factor(ytr)
#  }
  
#  return(list(mz=mz,logmz=logmz, mz.labels=mz.labels, 
#              xtr=xtr, ytr=iytr, sample.labels=data$sample.labels))
#}

ppc.subset.and.reshape<- function(data, user.parms){
  
 data.keep<-data

  xtr<-data$xtr
  mz<-data$mz
  logmz<-data$logmz
  ytr<-data$ytr
  
 mz.labels<-NULL

  sample.labels <- data$sample.labels
  patient.labels<- data$patient.labels
 
  if (!is.null(user.parms$mz.min) | !is.null(user.parms$mz.max)) {
    o1<- user.parms$mz.min
    o2<- user.parms$mz.max
    if(!is.null(o1)) {
      xtr<- xtr[mz>= o1,]
      logmz<- logmz[mz>= o1]
      mz<- mz[mz>= o1]
    }
    if (!is.null(o2)) {
      xtr<- xtr[mz<= o2,]
      logmz<- logmz[mz<= o2]
      mz<-    mz[mz<= o2]
    }
  }
  ## here I reshape the data into one column per patient. mz and logmz values are
  ## string out into one long vector, spaced apart

  if (!is.null(data$batch.labels)) {
    npatients<-length(unique(data$patient.labels))
    nbatches<-length(unique(data$batch.labels))
    p<-nrow(xtr)
    
    labelord<-order(data$batch.labels)
    xtr<-xtr[,labelord]
    if(!is.null(ytr)){ ytr<- ytr[labelord][1:npatients]}
    
    xtr3<-matrix(NA,nrow=nbatches*nrow(xtr),ncol=npatients)
    ii<-1
    for (i in 1:nbatches) {
      xtr3[ii:(ii+p-1),]<- xtr[, ((i-1)*npatients+1):(i*npatients)]
      ii<-(ii+p)
    }
    
    
    mz4<-NULL
    logmz<-NULL
# ROB fixed this bug (subtle error for pred test data)

    mz.labels<-NULL
#    fac<- 2*(max(mz)-min(mz))
#    faclog<- 2*(max(log(mz))-min(log(mz)))
    
fac<- 2*(user.parms$mz.max- user.parms$mz.min)
faclog<- 2*(log(user.parms$mz.max)- log(user.parms$mz.min+1))

    for(ii in 1:nbatches) {
      mz4<-c(mz4,mz+(ii-1)*fac)
      logmz<-c(logmz,log(mz)+(ii-1)*faclog)
      mz.labels<-c(mz.labels,paste(round(mz,1),"-",data$batches[ii],sep=""))
    }
    xtr<-xtr3
    mz<-mz4
    
    sample.labels<-data$patient.labels[labelord]
    sample.labels<-sample.labels[1:npatients]
    patient.labels<-sample.labels
   
  }
  
  iytr<-NULL
  if(!is.null(ytr)) {
    iytr<- as.factor(ytr)
  }
  
 data.keep$mz<-mz
 data.keep$logmz<-logmz
 data.keep$mz.labels<- mz.labels
 data.keep$xtr<-xtr
 data.keep$ytr<-iytr
 data.keep$sample.labels <- sample.labels
 data.keep$patient.labels <- patient.labels
 
 
  return( data.keep)
}


##
## Original by Rob
##
#ppc.subset.and.reshape.peakdata<- function(data, user.parms){


## subset and reshaping for data containing initial peak lists

                      
#peaklist<-data$peaklist
#logmz<-data$logmz
#mz<-data$mz
#ytr<-data$ytr

#mz.labels<-NULL
#n<-length(peaklist)

#if(!is.null(user.parms$mz.min) | !is.null(user.parms$mz.max)){
#	 o1<- user.parms$mz.min
#	 o2<- user.parms$mz.max
#	 if(!is.null(o1)){ 
#                      peaklist<-subset.peaklist(peaklist,log(o1),-1)
#                   logmz<- logmz[mz>= o1]
#                   mz<- mz[mz>= o1]
                  
#  	}
#  	if(!is.null(o2)){
#                      peaklist<-subset.peaklist(peaklist,log(o2),+1)
#                     logmz<- logmz[mz<= o2]
#                     mz<-    mz[mz<= o2]
#  	}
#}

## here I reshape the data into peaklist per patient. mz and logmz values are
## strong out into one long vector, spaced apart

#if(!is.null(data$batch.labels)){
#         nbatches<- length(unique(data$batch.labels))
#        npatients<-n/nbatches

#	p<-length(peaklist)

#	patient.ord<-order(data$patient.labels)
       
#	pk<-peaklist[patient.ord]
      
#         mz4<-NULL
#        logmz<-NULL
#        mz.labels<-NULL
#        fac<- 2*(max(mz)-min(mz))
#        faclog<- 2*(max(log(mz))-min(log(mz)))

#       pknew<-vector("list",npatients)
#        ii<-0
#	for(i in 1:npatients){
#           temp0<-NULL
#           for(j in 1:nbatches){
#             ii<-ii+1
#             temp<-pk[[ii]]
#             temp[,1]<-temp[,1]+faclog*(j-1)
#	     temp0<-rbind(temp0,temp)
#           }
#           pknew[[i]]<-temp0
#	}
	

#	mz.labels<-NULL
#        fac<- 2*(max(mz)-min(mz))
#        faclog<- 2*(max(log(mz))-min(log(mz)))

#	for(ii in 1:nbatches){
#	   mz4<-c(mz4,mz+(ii-1)*fac)
#	  logmz<-c(logmz,log(mz)+(ii-1)*faclog)
#	  mz.labels<-c(mz.labels,paste(round(mz,1),"-",data$batches[ii],sep=""))
#        }
#mz<-mz4
#if(!is.null(ytr)){  ytr<- matrix(ytr[patient.ord],ncol=nbatches,byrow=T)[,1]}
#sample.labels<- matrix(data$sample.labels[patient.ord],ncol=nbatches,byrow=T)[,1]
#peaklist<- pknew
#}

#iytr<-NULL
#if(!is.null(ytr)){ iytr<- as.factor(ytr)}

#return(list(mz=mz,logmz=logmz, mz.labels=mz.labels, peaklist=peaklist, ytr=iytr,
#sample.labels=sample.labels))
#}


"subset.peaklist"<- function(peaklist,r, dir){
       n<-length(peaklist)
       for(i in 1:n){
        if(dir== -1){
           o<-peaklist[[i]][,1]>r
              peaklist[[i]]<-peaklist[[i]][o,]
        }
     if(dir== +1){
           o<-peaklist[[i]][,1]<r
              peaklist[[i]]<-peaklist[[i]][o,]
        }
}
return(peaklist)
}
 ppc.subset.and.reshape.peakdata<-
function (data, user.parms) 
{
  data.keep <- data
    peaklist <- data$peaklist
    logmz <- data$logmz
    mz <- data$mz
    ytr <- data$ytr
    sample.labels <- data$sample.labels
    mz.labels <- NULL
    n <- length(peaklist)
    if (!is.null(user.parms$mz.min) | !is.null(user.parms$mz.max)) {
        o1 <- user.parms$mz.min
        o2 <- user.parms$mz.max
        if (!is.null(o1)) {
            peaklist <- subset.peaklist(peaklist, log(o1), -1)
            logmz <- logmz[mz >= o1]
            mz <- mz[mz >= o1]
        }
        if (!is.null(o2)) {
            peaklist <- subset.peaklist(peaklist, log(o2), +1)
            logmz <- logmz[mz <= o2]
            mz <- mz[mz <= o2]
        }
    }
    if (!is.null(data$batch.labels)) {
        nbatches <- length(data$batches)
        npatients <- n/nbatches
        p <- length(peaklist)
        patient.ord <- order(data$patient.labels)
        pk <- peaklist[patient.ord]
        mz4 <- NULL
        logmz <- NULL
        mz.labels <- NULL

# ROB fixed this bug (subtle error for test set pred!)

#        fac <- 2 * (max(mz) - min(mz))
#        faclog <- 2 * (max(log(mz)) - min(log(mz)))

        fac<- 2*(user.parms$mz.max- user.parms$mz.min)
        faclog<- 2*(log(user.parms$mz.max)- log(user.parms$mz.min+1))

        pknew <- vector("list", npatients)
        ii <- 0
        for (i in 1:npatients) {
            temp0 <- NULL
            for (j in 1:nbatches) {
                ii <- ii + 1
                temp <- pk[[ii]]
                temp[, 1] <- temp[, 1] + faclog * (j - 1)
                temp0 <- rbind(temp0, temp)
            }
            pknew[[i]] <- temp0
        }
        mz.labels <- NULL
        fac <- 2 * (max(mz) - min(mz))
        faclog <- 2 * (max(log(mz)) - min(log(mz)))
        for (ii in 1:nbatches) {
            mz4 <- c(mz4, mz + (ii - 1) * fac)
            logmz <- c(logmz, log(mz) + (ii - 1) * faclog)
            mz.labels <- c(mz.labels, paste(round(mz, 1), "-", 
                data$batches[ii], sep = ""))
        }
        mz <- mz4
        if (!is.null(ytr)) {
            ytr <- matrix(ytr[patient.ord], ncol = nbatches, 
                byrow = T)[, 1]
        }
        sample.labels <- matrix(data$sample.labels[patient.ord], 
            ncol = nbatches, byrow = T)[, 1]
        peaklist <- pknew
    }
    iytr <- NULL
    if (!is.null(ytr)) {
        iytr <- as.factor(ytr)
    }
    data.keep$mz <- mz
    data.keep$logmz <- logmz
    data.keep$mz.labels <- mz.labels
    data.keep$peaklist <- peaklist
    data.keep$ytr <- iytr
    data.keep$sample.labels <- sample.labels
    return(data.keep)
}

##
## Functions for the Excel interface
##
##
## Read all data files in a folder and return mz and the matrix of spectra
##

##
## A list that is used for setting ppc.options
##
ppc.options <- list(debug=FALSE, #whether to turn on debugging or not
                    err.file=ifelse(.Platform$OS.type=="windows", "C:/ppctrace.txt", "ppctrace.txt"),
                    reserved.class.label="Unspecified")

ppc.constants <- list(data.types = list(raw="RAW", peak="PEAK"))


##
## Our error handler
##
ppc.xl.error.trace <- function() {
  err.message <- geterrmessage()
  sink(ppc.options$err.file)
  print(err.message)
  traceback()
  sink()
  winDialog(type="ok", message=err.message)
}

##
## Upon loading, if we are in a windows environment, we use the windows
## dialog mechanism to display errors. Useful for debugging COM apps
##
.onLoad <- function(lib, pkg) {
  if ( .Platform$OS.type == "windows") {
#    options(error=function() winDialog(type="ok", message=geterrmessage()))
    options(error=ppc.xl.error.trace)
  }

}

##
## Upon unload, we set things back the way they were...
##
.onUnload <- function(libpath){
  if ( .Platform$OS.type == "windows") {
    options(error=NULL)
  }
}
       
##
## A function to get the names of subdirs in a given directory.
## Note: full pathnames are returned.
##
ppc.xl.get.names.of.subdirs <- function(directory) {
  f.info <- file.info(dir(directory, full.names=TRUE))
  return(dimnames(f.info)[[1]][f.info[, "isdir"] == TRUE])
}

##
## A function to get the names of files in a given directory.
## Note: full pathnames are returned.
##
ppc.xl.get.names.of.files <- function(directory) {
  f.info <- file.info(dir(directory, full.names=TRUE))
  return(dimnames(f.info)[[1]][f.info[, "isdir"] == FALSE])
}

##
## NOTE to memory-leaking myself:
##       All the functions below expect the class data directory
##       that is, the directory that is the root of all class specific data
##       We have to call the functions for each class. The plurals expect the entire list.
##

##
## A function to get detect if we have batches or not.
##
ppc.xl.detect.batches <- function(classDataDirectory) {
  x <- ppc.xl.get.names.of.subdirs(classDataDirectory)
  return(length(x) >= 2)  # Need at least two batches!
}

##
## A function to get detect the batch labels
## Note: only the batch labels (not full pathname) are returned
##       Also, batches are assumed to be present.
##       Batch labels are merely base names of subdirectories 
##
ppc.xl.compute.batch.labels <- function(classDataDirectories) {
  x <- unlist(sapply(classDataDirectories,
                     function(x) ppc.xl.get.names.of.files(ppc.xl.get.names.of.subdirs(x))))
  return (as.vector(sapply(x, function(x) basename(dirname(x)))))
}

##
## A function to get the list of all data files 
## Note: Full path names are returned
##
ppc.xl.get.all.data.file.names <- function(classDataDirectories, batches.present=FALSE) {
  if (batches.present) {
    x <- unlist(sapply(classDataDirectories,
                       function(x) ppc.xl.get.names.of.files(ppc.xl.get.names.of.subdirs(x))))
  } else {
    x <- unlist(sapply(classDataDirectories,
                       function(x) ppc.xl.get.names.of.files(x)))
  }
  return(as.vector(x))
}

##
## A function to compute sample labels
## Note: if batches are present, sample labels are <batch_name>-<file-name-without-suffix>
##       otherwise just <file-name-without-suffix>
##
ppc.xl.compute.sample.labels <- function(classDataDirectories, batches.present=FALSE) {
  if (batches.present) {
    x <- unlist(sapply(classDataDirectories,
                       function(x)
                       basename(ppc.xl.get.names.of.files(ppc.xl.get.names.of.subdirs(x)))))
    x <- paste(ppc.xl.compute.batch.labels(classDataDirectories), x, sep="-")
  } else {
    x <- unlist(sapply(classDataDirectories,
                       function(x)
                       basename(ppc.xl.get.names.of.files(x))))
  }
  return(as.vector(sapply(x, ppc.xl.remove.suffix)))
}

##
## A function to compute sample sizes
## Note: Input is the vector of data folders
##
ppc.xl.compute.sample.sizes <- function(classDataDirectories, batches.present=FALSE) {
  result <- sapply(classDataDirectories,
                   function(x)
                   ifelse(batches.present,
                          length(ppc.xl.get.names.of.files(ppc.xl.get.names.of.subdirs(x))),
                          length(ppc.xl.get.names.of.files(x))))
  return (as.vector(result))
}


##
## A function to compute patient labels
## Note: patient labels are merely basenames of files without the extension
##
ppc.xl.compute.patient.labels <- function(classDataDirectories, batches.present=FALSE) {
  if (batches.present) {
    x <- unlist(sapply(classDataDirectories,
                       function(x)
                       basename(ppc.xl.get.names.of.files(ppc.xl.get.names.of.subdirs(x)))))
  } else {
    x <- unlist(sapply(classDataDirectories,
                       function(x)
                       basename(ppc.xl.get.names.of.files(x))))
  }
  return(as.vector(sapply(x, ppc.remove.suffix)))
}


##
## A function to get the names of files in a given directory for a given batch
## Note: full pathnames are returned
##
ppc.xl.get.files.in.batch <- function(classDataDirectory, batchName) {
  f.info <- file.info(dir(file.path(classDataDirectory, batchName), full.names=TRUE))
  return(dimnames(f.info)[[1]][f.info[, "isdir"] == FALSE])
}



##
## Read data from a single file. 
##
## 
ppc.xl.read.data.file  <- function(file.name, mz=NULL) {
  if (ppc.options$debug) {
    print(paste("Reading file", file.name))
  }

  file.contents <- read.table(file.name, sep=",")

  if (!is.null(mz)) {
    y <- approx(file.contents[, 1], file.contents[, 2], xout = mz)$y
  } else {
    mz <- file.contents[, 1]
    y <- file.contents[, 2]
  }
  return(list(mz=mz, y = y))
}

##
## Build raw data set
##
ppc.xl.build.data  <- function(raw.file.data,
                               class.labels,
                               sample.labels,
                               batch.labels,
                               patient.labels,
                               batches.exist = FALSE,
                               data.type=ppc.constants$data.types$raw,
                               class.levels = NULL) {
  if (is.null(class.levels)) {
    ytr <- factor(class.labels)
  } else {
    ytr <- factor(class.labels, levels=class.levels)
  }

  batches <- NULL
  if (batches.exist) {
    batches <- unique(batch.labels)
  }
  
  if (data.type == ppc.constants$data.types$raw) {
    return(list(mz=raw.file.data$mz,
                logmz = log(raw.file.data$mz),
                xtr = raw.file.data$xtr,
                ytr = ytr,
                peaklist = raw.file.data$peaklist, ## will be NULL!
                sample.labels = sample.labels,
                batch.labels = batch.labels,
                batches = batches,
                patient.labels = patient.labels,
                data.type = data.type))
  } else {
    return(list(mz=raw.file.data$mz,
                logmz = log(raw.file.data$mz),
                xtr = NULL,
                ytr = ytr,
                peaklist = raw.file.data$peaklist,
                sample.labels = sample.labels,
                batch.labels = batch.labels,
                batches = batches,
                patient.labels = patient.labels,
                data.type = data.type))
  }
}

##
## Get the current user parameters
##
ppc.xl.get.user.parameters  <- function() {
  return(ppc.xl.current.user.parameters)
}

##
## Set the current user parameters
##
ppc.xl.set.user.parameters  <- function(x) {
  ppc.xl.current.user.parameters  <- x
}

##
## Return the default set of parameters
##
ppc.xl.get.default.parameters  <- function() {
  return (list(stn = 1,
               minht = 0,
               span = 201,
               smoothing.span = 0.05,
               nsplits = 10,
               fix.at.one = FALSE,
               peak.gap = 0.005,
               recluster = FALSE,
               mz.min = 15,
               mz.max = 1500,
               nperms=20))
}

##
## Construct the current set of parameters
##
ppc.xl.make.parameter.set  <- function(stn=1, minht=0, span=201,
                                       smoothing.span=0.05,
                                       nsplits=10,
                                       fix.at.one=FALSE,
                                       peak.gap=0.005,
                                       recluster=FALSE,
                                       mz.min=15,
                                       mz.max=1500,
                                       nperms=20) {
  return (list(stn = stn,
               minht = minht,
               span = span,
               smoothing.span = smoothing.span,
               nsplits = nsplits,
               fix.at.one = fix.at.one,
               peak.gap = peak.gap,
               recluster = recluster,
               mz.min = mz.min,
               mz.max = mz.max,
               nperms=nperms))
}


##
## Compute the training errors for the dataset
ppc.xl.compute.training.errors  <- function(fit, data) {
  foo  <- ppc.error(data$ytr, fit$yhat)
  threshold = fit$threshold
  n  <- length(data$ytr)
  
  return (list(x = fit$threshold,
               y = foo$err/n,
               y.ytop = fit$numsites,
               x.label = "Threshold",
               y.label = "Training Error"))
}

##
## Compute the test errors for the dataset
## TODO Need to fix.
ppc.xl.compute.test.errors  <- function(fit, data) {
  err.cnts  <- ppc.error(data$ytr, fit$yhat)
  threshold <- fit$threshold
  n  <- length(data$ytr)
  
  return (list(x = fit$threshold,
               y = err.cnts$err/n,
               y.ytop = fit$numsites,
               x.label = "Threshold",
               y.label = "Test Error"))
}


##
## Rob likes to predict a whole set of things in one shot, but there is nothing that says
## that a user might not want to predict for a particular value of a threshold.
## SO I have to rewrite the predict.ppc function to do things just for a single threshold.
##
ppc.xl.predict <- function(centroids.fit, split.fit, logmz, peaklist.te, threshold=NULL,
                           metric=c("binomial","euclidean", "absolute"),
                           summ=c("mean","median")) {
  
  ## test set prediction for PPC method
  ## makes  predictions  for test set list of peaks in peaklist.te
  ##  
  ## "metric" is the metric used in the nearest centroid rule; "summ" is
  ## summary used  to combine distances over sites
  ## n.threshold is the number of shrinkage thresholds used
  
  ## returns yhat, plus ind0= indicator matrix  of whether peak was found,
  ##           ind=indicator matrix  of whether peak > cutpoint was found
  ##           ht=matrix of peak heights
  
  ## "numsites" is the number of sites present at each threshold; "sites" is a list
  ## containing the ##s of these sites
  ##
  ## dis is the distance of each test profile to each class centroid
  
  ##
  metric <- match.arg(metric)
  summ <- match.arg(summ)
  ##  print(metric)
  ##  print(summ)
  
  n <- length(peaklist.te)
  n.class <- split.fit$n.class
  K <- ncol(split.fit$prhat)
  p <-length(split.fit$cuthat)
  yhat <- matrix(NA, nrow = p, ncol = n)
  yhatt <- vector(mode="numeric", length = n.class)
  ind0 <- yhat
  ht <- yhat
  dis <- array(NA,c(n, K))

  for(j in 1:n) {
    if (ppc.options$debug) cat(j)
    aaa <- ppc.predict.peaks1(centroids.fit, logmz,  peaklist= 
                              peaklist.te[[j]])
    ht.hat <-  aaa$ht * aaa$ind
    ind0[,j] <- 1*aaa$ind
    ht[,j] <- ht.hat
    for(i in 1:p) {
      yhat[i, j] <- 1 * (ht.hat[i] > split.fit$cutpoints[i, split.fit$cuthat[i]])
    }
  }        
  
  prmean <- apply(split.fit$prhat,1,mean)
  delta <- split.fit$prhat-matrix(prmean,nrow=nrow(split.fit$prhat),ncol=K)
  
  dd<- matrix(NA,nrow=n,ncol=K)
  
  ##  if(is.null(threshold)){threshold <- seq(0, max(abs(delta)), length = n.threshold)}
  
  soft.thresh <- function(x, tt) {
    sign(x) * (abs(x) - tt) * (abs(x) > tt)
  }
  
  delta2 <- soft.thresh(delta, threshold)
  pr <- prmean + delta2
  sumabs<-  apply(abs(delta2),1,sum)
  pos <-  sumabs!=0
  numsites <- sum(pos)
  sites <- (1:nrow(pr))[pos]
  
  sites <- sites[order( - abs(sumabs[pos]))]
  if(metric == "euclidean") {
    for(k in 1:K){
      temp<- (yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))^2
      dd[,k] <- apply(temp,2,summ)
    }
  }
  if(metric == "absolute") {
    for(k in 1:K){
      temp<- abs(yhat - matrix(pr[,k], ncol = n, nrow = nrow(yhat)))
      dd[,k] <- apply(temp,2,summ)
    }
  }
  if(metric == "binomial") {
    ## compute Bayes estimate of pr, to keep it away from 0 or 1
    ntemp <- matrix(n.class,nrow=nrow(pr),ncol=length(n.class), byrow=T)
    prbayes<- (pr*ntemp+1)/(ntemp+2)
    for(k in 1:K){
      prmat<-matrix(prbayes[,k], ncol = n, nrow = nrow(yhat))
      temp<- -1*(yhat*log(prmat)+(1-yhat)*log(1-prmat))
      dd[,k] <- apply(temp,2,summ)
    }
  }
  
  yhatt <- apply(dd,1,which.is.min)
  dis <- dd

  prob <- exp(-dis)
  prob <- prob/apply(prob,1,sum)
  
  return(list(yhat = yhatt, threshold = threshold, numsites = numsites,
              sites = sites, ind = yhat, ind0=ind0, prob=prob, ht=ht))
  
}


##
## Compute the confusion matrix for the dataset
ppc.xl.compute.confusion.matrix  <- function(centroids.fit, split.fit,
                                             data, peaklist.te, threshold=NULL,
                                             metric=c("binomial","euclidean", "absolute"),
                                             summ=c("mean","median")) {
  
  ##i has to be determined!
  metric  <- match.arg(metric)
  summ  <- match.arg(summ)
  pred.obj  <- ppc.xl.predict(centroids.fit, split.fit,
                              data$logmz, peaklist.te=peaklist.te, threshold=threshold,
                              metric=metric, summ=summ)
  
  mat  <- table(data$ytr, pred.obj$yhat)
  return(list(confusion.matrix=mat))
}

##
## massage arrays and matrices to handle NAs in Excel
##
ppc.xl.transform.matrix  <- function(x) {
  w  <- x
  m  <- apply(w, 2, function(x) {ifelse(is.na(x),1,0)})
  w  <- apply(w, 2, function(x) {ifelse(is.na(x),0,x)})
  return(list(matrix=w, missing=m))
}

##
## Compute the ppc cv errors
##
ppc.xl.compute.cv.errors  <- function(trained.obj, cv.obj) {
  n <- nrow(cv.obj$yhat)
  y <- cv.obj$y
  nc <- length(table(y))
  nfolds <- length(cv.obj$folds)
  err2 <- matrix(NA, nrow = length(unique(y)), ncol = length(cv.obj$threshold))
  for(i in 1:(length(cv.obj$threshold))) {
    s <- cv.obj$confusion[[i]]
    diag(s) <- 0
    err2[, i] <- apply(s, 1, sum)/table(y)
  }
  
  return(list(x=trained.obj$threshold, y=cv.obj$err,
              x.label="Threshold", y.label="Error",
              y.se=cv.obj$se, numsites=cv.obj$numsites,
              cv.err=t(err2), cv.legend= dimnames(table(y))[[1]]))
}

##
## Find the nearest value to a list of nearest value
##
ppc.xl.find.index.of.nearest.value  <- function(target.vec, val) {
  return(which(order(abs(target.vec - val)) == 1))
}


##
## Setup for cv
##
ppc.xl.cv.setup <- function(ppc.fit, data,   user.parms) {
  
  ## K-fold cross-validation for ppc method
  ## takes results of make.centroids.list, find.splits, and predict.ppc
  ## and computes cross-validated predictions and error estimate.

  peaklist.tr <- data$peaklist
  ytr <- data$ytr
  logmz <- data$logmz
  
  peak.gap <- user.parms$peak.gap
  nsplits <- user.parms$nsplits
  fix.at.one <- user.parms$fix.at.one
  recluster <- user.parms$recluster
  
  folds <- balanced.folds(ytr)
  
  n.class<-table(ytr)
  n.threshold<- length(ppc.fit$threshold)
  
  n <- length(peaklist.tr)
  yhatcv <- array(NA,c(n,n.threshold))
  probcv <- array(1, c(n, length(n.class), n.threshold))
  return(list(numsites = ppc.fit$numsites,
              peaklist.tr = peaklist.tr,
              ytr = ytr,
              logmz = logmz,
              peak.gap = peak.gap,
              nsplits = nsplits,
              fix.at.one = fix.at.one,
              recluster = recluster,
              folds = folds,
              threshold = ppc.fit$threshold,
              n.threshold = n.threshold,
              n = n,
              yhatcv = yhatcv,
              probcv = probcv))
}

##
## DO fold ii
##
ppc.xl.cv.do.fold <- function(cv.state, ii, user.parms) {
  gg <- cv.state$folds[[ii]]
  if (ppc.options$debug) print(c("fold=",ii),fill=T)
  pk <- cv.state$peaklist.tr[-gg]
  data.temp <- list(ytr=cv.state$ytr[-gg], logmz=cv.state$logmz, peaklist=pk)
  astar <- ppc.make.centroid.list(data.temp,  user.parms)
  junk0 <- ppc.predict.peaks(astar,  data.temp)
  aa <- ppc.find.splits(astar, junk0 ,data.temp,user.parms)
  ress <- ppc.predict(astar,aa,cv.state$logmz,peaklist.te=cv.state$peaklist.tr[gg],
                      threshold=cv.state$threshold, metric="euclidean")
  return(list(gg=gg, yhat=ress$yhat, prob=ress$prob))
  ##  yhatcv[gg,] <- ress$yhat
  ##  cv.state$probcv[gg,,] <- ress$prob
}

##
## Compute the cv probabilities
##
ppc.xl.compute.cvprobs  <- function(cv.fit, data, threshold) {
  ii <- (1:length(cv.fit$threshold))[cv.fit$threshold > threshold]
  ii <- ii[1]
  ss <- data$sample.labels
  pp <- cv.fit$prob[,  , ii]
  
  y <- cv.fit$y
  
  o <- order(y)
  y <- y[o]
  if(!is.null(ss)) {
    ss <- ss[o]
  }
  ppp <- pp[o,  ]
  n <- nrow(ppp)
  nc <- length(unique(y))
  
  return (list(x = 1:n,
               y = ppp,
               x.label = "Sample",
               y.label = "CV Probabilities",
               y.names = levels(factor(y)),
               y.lines = cumsum(table(cv.fit$y)),
               x.dummy = vector(length=2, mode="numeric"),
               y.dummy = vector(length=2, mode="numeric"),
               x.names = ss))
}


##
## Compute the prediction probabilities for test set.
##
  
ppc.xl.compute.test.probs  <- function(centroid.fit, split.fit, training.class.names, data, peaklist,
                                       threshold) {
  ppc.fit1 <-  ppc.predict1(centroid.fit, split.fit, data$logmz, peaklist, threshold)
  predicted.y <- sapply(ppc.fit1$yhat, function(x){ training.class.names[x] })
  predicted.probs <- ppc.fit1$prob
  sample.labels <- data$sample.labels

  order.classes  <- order(predicted.y)
  actual.classes <- as.character(data$ytr[order.classes])
  actual.classes[is.na(actual.classes)] <- ppc.options$reserved.class.label
  pp <- apply(predicted.probs, 2, function(x) x[order.classes])
  ny  <- predicted.y[order.classes]
  n  <- length(ny)
  sample.labels  <- data$sample.labels[order.classes]
  actual.class.names <- levels(factor(actual.classes))

  return (list(x = 1:n,
               y = pp,
               x.label = "Sample",
               y.label = "Predicted Test Probabilities",
               y.names = training.class.names,
               y.lines = cumsum(table(actual.classes)) + 0.5,
               x.dummy = vector(length=2, mode="numeric"),
               y.dummy = vector(length=2, mode="numeric"),
               panel.names = actual.class.names,
               x.names = sample.labels))
}

##
## Compute the quantities for the ``Show Prediction'' button in PPC (Excel)
##

ppc.xl.compute.test.prediction  <- function(centroid.fit, split.fit, training.class.names, data, peaklist,
                                       threshold) {
  ppc.fit1 <-  ppc.predict1(centroid.fit, split.fit, data$logmz, peaklist, threshold)
  predicted.y <- sapply(ppc.fit1$yhat, function(x){ training.class.names[x] })
  predicted.probs <- ppc.fit1$prob
  colnames(predicted.probs) <- training.class.names
  actual.classes <- as.character(data$ytr)
  data.has.missing.class.labels <- any(is.na(actual.classes))
  actual.classes[is.na(actual.classes)] <- ppc.options$reserved.class.label
  
  if (data.has.missing.class.labels) {
    confusion.matrix <- NULL
  } else {
    confusion.matrix <- table(actual.classes, predicted.y)
  }
  return (list(actual.class.labels = actual.classes,
               predicted.class.labels = predicted.y,
               sample.labels = data$sample.labels,
               patient.labels = data$patient.labels,
               confusion.matrix = confusion.matrix,
               predicted.probs = t(predicted.probs)))
}



##
## Transform class labels by replacing NAs with reserved label
##
ppc.xl.transform.class.labels <- function(class.labels, reservedLabel) {
  w <- class.labels
  w[is.na(class.labels)] <- reservedLabel
  return (w)
}

##
## Plotting of histograms
##
ppc.xl.plothist <-  function(peak.fit, ppc.fit, centroid.fit, split.fit, data, nsites=NULL) {

  require(class)
  sitelist <- ppc.fit$sites[[1]]
  
  if(!is.null(nsites)){
    sitelist<-sitelist[1:nsites]
  }
  
  ytr<- data$ytr
  codesy <- levels(ytr)
  
  nc<- dim(split.fit$pr)[3]
  nrows <- trunc(20/nc)
  win.metafile()
  
  par(mfcol=c(nrows,6))
  
  par(mar=c(0,1,0,1))
  par(cex=.6)
  ht <- peak.fit$ind*peak.fit$ht
  
  
  for(i in sitelist){
    
    if (ppc.options$debug) cat(i,fill=T)
    xmax <- max(ht[i,])
    br <- hist(ht[i,],col=3,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=F, main="")$br
    
    par(mar=c(0,1,1,1))
    junk <- hist(ht[i,ytr==codesy[1]],col=3,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=F
                 ,breaks=br)
    hist(ht[i,ytr==codesy[1]],col=5,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br, main="")
    box(col=gray)
    wid <- (junk$br[2]-junk$br[1])
    v <- split.fit$cutpoints[i,split.fit$cuthat[i]]+wid
    v0<-v
    if(v==0){v<-v+wid}
    
    h <- max(junk$counts)
    
    abline(v=v,lty=1,col=2)
    
    for(j in (1:ncol(split.fit$prclose))[split.fit$prclose[i,]]){
      points(split.fit$cutpoints[i,j]+wid,h/2,col=2,pch="x")
    }
    
    o<-knn1(data$logmz, centroid.fit$cent[i,1], 1:length(data$logmz))
    lab<- data$mz.labels[o]
    text(xmax*.5,.9*h,labels=lab,cex=.8)
    text(v,3*h/4,labels=as.character(round(v0,2)),col=6)
    
    
    pr <- sum(ht[i,ytr==codesy[1]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[1])
    text(xmax,h/2,labels=as.character(round(pr,2)),col=2,cex=.6)
    
    if(nc>2){
      for(ii in 2:(nc-1)){
        par(mar=c(1,1,1,1))
        junk2 <- hist(ht[i,ytr==codesy[ii]],col=4,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=
                      F,breaks=br)
        hist(ht[i,ytr==codesy[ii]],col=ii,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br, main="") 
        
        box(col=gray)
        abline(v=v,lty=1,col=2)
        h2 <- max(junk2$counts)
        pr <- sum(ht[i,ytr==codesy[ii]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[ii])
        text(xmax,h2/2,labels=as.character(round(pr,2)),col=2,cex=.6)
      }}
    
    par(mar=c(2,1,0,1))
    junk2 <- hist(ht[i,ytr==codesy[nc]],col=4,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,plot=
                  F,breaks=br)
    hist(ht[i,ytr==codesy[nc]],col=nc,xlab="",nclass=10,xlim=c(0,xmax*1.15),axes=F,breaks=br,main="") 
    
    box(col=gray)
    abline(v=v,lty=1,col=2)
    h2 <- max(junk2$counts)
    pr <- sum(ht[i,ytr==codesy[nc]]>split.fit$cutpoints[i,split.fit$cuthat[i]])/sum(ytr==codesy[nc])
    text(xmax,h2/2,labels=as.character(round(pr,2)),col=2,cex=.6)
    
  }
  dev.off()
  return()
}
  
##
## Read a single peaks file and return list of logmz and peak
##
ppc.xl.read.peak.file <- function(filename) {
  ## note: we ignore anything wiht mz <=0 
  pat<-scan(filename,what="", comment.char="#")
  pat<-matrix(pat, ncol=7, byrow=T) [,(1:2)]
  pat<-matrix(as.numeric(as.character(pat)),ncol=2)
  pat <- pat[pat[, 1] > 0, ]
  # rob added log in the next line
  return(cbind(log(pat[,1]),pat[,2]))
}

##
## Read peak data from given files.
##
ppc.xl.read.peak.data <- function(files, batches.exist) {
  ##  if (! batches.exist) {
  n <- length(files)
  peaklist <- vector("list",n)
  mz <- NULL
  for(i  in 1:n){
    w <- ppc.xl.read.peak.file(files[i])
    peaklist[[i]] <- w
    # rob added the exp in the next line
    mz <- c(mz, exp(w[, 1]))
  }
  return(list(peaklist=peaklist, mz=mz))
  ##}
}

##
## Subset and reshape data according to its type
##
ppc.xl.subset.and.reshape <- function(data, user.params) {
  if (data$data.type==ppc.constants$data.types$raw) {
    return(ppc.subset.and.reshape(data, user.params))
  } else {
    return(ppc.subset.and.reshape.peakdata(data, user.params))
  }
}

##
## Depending on the dataset, augment the data with peaks if necessary
##
ppc.xl.augment.with.peaks <- function(data, user.parms) {
  if (data$data.type==ppc.constants$data.types$raw) {
    ## add peaklist to data object
    data$peaklist<- ppc.make.peaklist(data, user.parms)
  }
  return(data)
}


ppc.xl.plotfdr <- function(data, centroid.fit, peak.fit, split.fit, training.fit, user.parms) {
  fdrfit <- ppc.fdr(data, centroid.fit, peak.fit, split.fit, training.fit, user.parms)
  return(list(x=fdrfit$results[,"npeaks"],
              y=fdrfit$results[,"fdr"],
              x.label="No. of significant peaks",
              y.label="False discovery rate"))
}

ppc.xl.remove.suffix <- function(x){
  w <- unlist(strsplit(x, "\\."))
  m <- length(w)
  return(ifelse(m==1, w, paste(w[1:(m-1)], sep="", collapse=".")))
}

