.packageName <- "superpc"
superpc.cv <- function (fit, data,  n.threshold = 20, 
                        n.fold = 10, folds=NULL,  n.components=3,
                        min.features=5, max.features=nrow(data$x), compute.fullcv=FALSE)
  
 # cross-validation for supervised PCs;
 # returns both  preval cv and fullcv (if requested)  

  
{
  this.call <- match.call()
  type <- fit$type

  if(n.components>5){ cat("Max # of components is 5",fill=TRUE)}

  n.components <- min(5, n.components)

  mean.na<- function(x){mean(x[!is.na(x)])}

  n <- ncol(data$x)
  cur.tt <- fit$feature.scores

  lower <- quantile(abs(cur.tt), 1 - (max.features/nrow(data$x)))
  upper <- quantile(abs(cur.tt), 1 - (min.features/nrow(data$x)))


  if(is.null(folds)){
    folds<-vector("list",n.fold)
    breaks <- round(seq(from = 1, to = (n + 1), length = (n.fold +
                          1)))
    cv.order <- sample(1:n)
    for(j in 1:n.fold){
      folds[[j]]<-cv.order[(breaks[j]):(breaks[j + 1] - 1)]
    }
  }

  featurescores.folds <- matrix(nrow=nrow(data$x), ncol= n.fold)

  th <- seq(from = lower, to = upper, length = n.threshold)
  nonzero <- rep(0, n.threshold)
  out <- array(NA, c(n.components, n.threshold, n.fold))
  out.preval<-matrix(NA,nrow=n.components, ncol=n.threshold)

  cur2<-array(NA,c(n,n.components,n.threshold))

# note, unlike in superpc.predict, we do not flip the signs of the latent
#factors. I don;t think this will cause a problem!


  for (j in 1:n.fold) {
    cat("",fill=TRUE)
    cat(c("fold=",j),fill=TRUE)
    data.temp=list(x=data$x[,-folds[[j]]], y=data$y[-folds[[j]]], status=data$status[-folds[[j]]])
    cur.tt <- superpc.train(data.temp, type=type)$feature.scores
    featurescores.folds[,j]<- cur.tt
    for (i in 1:n.threshold) {
      cat(i)
      cur.features <- (abs(cur.tt) > th[i])
      if(sum(cur.features)>1){
        
        nonzero[i] <-  nonzero[i]+sum(cur.features)/n.fold
        
        
        cur.svd <- mysvd(data$x[cur.features, -folds[[j]]],n.components=n.components)
        

        cur.v.all <- scale(t(data$x[cur.features, folds[[j]], drop=FALSE]) %*% cur.svd$u, center=FALSE,scale=cur.svd$d)

        n.components.eff<- min(sum(cur.features),n.components)
        cur.v  <- cur.v.all[,1:n.components.eff]
        cur2[folds[[j]],1:n.components.eff, i]<-cur.v

        if(compute.fullcv){
          for (k in 1:ncol(cur.v)) {
            if(type=="survival"){
              junk <- coxph(Surv(data$y[folds[[j]]], data$status[folds[[j]]]) ~cur.v[, 1:k])$loglik
              out[k,i, j]<-2*(junk[2]-junk[1]) 
            }
            else{
              junk<-summary(lm(data$y[folds[[j]]]~cur.v[, 1:k]))
              out[k,i,j]<- junk$fstat[1]
            }

          }
        }
      }
    }
  }
  cat("\n")

  out<- apply(out,c(1,2),mean.na)

  for(i in 1:n.threshold){
    for(j in 1:n.components){
      if(type=="survival"){
        require(survival)
        junk<-  coxph(Surv(data$y, data$status) ~cur2[, 1:j, i])$loglik
        out.preval[j,i]<- 2*(junk[2]-junk[1])
      }
      else{
        junk<-summary(lm(data$y~cur2[, 1:j, i]))
        out.preval[j,i]<- junk$fstat[1]
      }

    }}

  junk <- list(threshold = th, nonzero=nonzero, scor.preval=out.preval, scor=out,
               folds=folds, featurescores.folds=featurescores.folds,  v.preval=cur2,  type=type, call = this.call)
  class(junk) <- "superpc.cv"
  return(junk)
}


superpc.fit.to.outcome<- function(fit, data.test,score, print=TRUE){


type=fit$type

 if(type=="survival"){
   require(survival)
   result<-coxph(Surv(data.test$y, data.test$status)~score)
}

 else{
   result<-lm(data.test$y~score)
}


if(print){print(summary(result))}

return(result)
}

superpc.listfeatures<- function(data, train.obj, fitred, component.number, shrinkage){

if( shrinkage < min(fitred$shrinkages) |  shrinkage > max(fitred$shrinkages)){
    stop("Error: shrinkage value out of range")
}

temp<- abs(shrinkage- fitred$shrinkages)

ii<-(1:length(temp))[temp==min(temp)]

featurenames.short<- substring(data$featurenames,1,40)

oo=fitred$feature.list[[ii]][[component.number]]

res<-cbind(round(fitred$import[oo,component.number],3), round(train.obj$feature.scores[oo],3), featurenames.short[oo])



o<-order(-abs(fitred$import[oo,component.number]))
res<-res[o,]
dimnames(res)<-list(NULL,c("Importance-score", "Raw-score", "Name"))
return(res)
}
superpc.lrtest.curv <- function (object, data, newdata, n.components=1, threshold=NULL,  n.threshold=20)
{

  this.call <- match.call()

                                        # compute lrtest statistics based on fit "object", training data "data",
                                        # and test  data "newdata",
                                        # over a set of threshold values

  type=object$type

  if(!is.null(threshold)) {n.threshold=length(threshold)}
  if(is.null(threshold)){
    second.biggest<- -sort(-abs(object$feature.scores))[2]
    threshold<- seq(0,second.biggest, length=n.threshold)
  }

  n.pc <- n.components
  lrtest<-rep(NA, n.threshold)
  num.features<-rep(NA, n.threshold)


  cat("",fill=TRUE)
  for(ii in 1:n.threshold){
    cat(ii)

    object.temp<- superpc.predict(object, data, newdata,threshold=threshold[ii], n.components=n.pc)
    
    num.features[ii]<- sum(object.temp$which.features)
    
    v.pred<-object.temp$v.pred
    
    
    if(type=="survival"){
      require(survival)
      junk<- coxph(Surv(newdata$y, newdata$status) ~v.pred)$loglik
      lrtest[ii]<-2*(junk[2]-junk[1])
    }
    else{junk<- summary(lm(newdata$y~v.pred))
         lrtest[ii]<-junk$fstat[1]
       }

  }
  cat("",fill=TRUE)

  return(list(lrtest=lrtest,threshold=threshold,num.features=num.features, type=type, call=this.call))
}
cor.func<- 

function (x, y, fudge = median(sd)) 
{
    n <- length(y)
    xbar <- x %*% rep(1/n, n)
    sxx <- ((x - as.vector(xbar))^2) %*% rep(1, n)
    sxy <- (x - as.vector(xbar)) %*% (y - mean(y))
    syy <- sum((y - mean(y))^2)
    numer <- sxy/sxx
    sd <- sqrt((syy/sxx - numer^2)/(n - 2))
    tt <- numer/(sd + fudge)
    return(list(tt = tt, numer = numer, sd = sd))
}

coxfunc <- 
function(x, y, status, fudge = median(sd))
{
        junk <- coxscor(x, y, status)
        scor<-junk$scor
        sd <- sqrt(coxvar(x, y, status, coxstuff.obj=junk$coxstuff.obj))
        tt <- scor/(sd + fudge)
        return(list(tt = tt, numer = scor, sd = sd))
}


coxscor <- 
function(x, y, ic, offset = rep(0., length(y)))
{
        # computes cox scor function for rows of nx by n matrix  x
 # first put everything in time order
        n <- length(y)
        nx <- nrow(x)
        yy <- y + (ic == 0.) * (1e-05)
        otag <- order(yy)
        y <- y[otag]
        ic <- ic[otag]
        x <- x[, otag, drop = F]
        #compute  unique failure times, d=# of deaths at each failure time, 
        #dd= expanded version of d to length n, s=sum of covariates at each
        # failure time, nn=#obs in each risk set, nno=sum(exp(offset)) at each failure time
        offset <- offset[otag]
        a <- coxstuff(x, y, ic, offset = offset)
        nf <- a$nf
        fail.times <- a$fail.times
        s <- a$s
        d <- a$d
        dd <- a$dd
        nn <- a$nn
        nno <- a$nno
        w <- rep(0., nx)
        for(i in (1.:nf)) {
                w <- w + s[, i]
                oo<- (1.:n)[y >= fail.times[i]]
                r<-rowSums(x[, oo, drop = F] * exp(offset[oo]))
                w<- w - (d[i]/nno[i])*r 
        }
        return(list(scor = w, coxstuff.obj = a))
}



 coxvar <- 
function(x, y, ic, offset = rep(0., length(y)), coxstuff.obj = NULL)
{
        # computes information elements (var) for cox
        # x is nx by n matrix of expression  values
        nx <- nrow(x)
        n <- length(y)
        yy <- y + (ic == 0.) * (1e-06)
        otag <- order(yy)
        y <- y[otag]
        ic <- ic[otag]
        x <- x[, otag, drop = F]
        offset <- offset[otag]
        if(is.null(coxstuff.obj)) {
                coxstuff.obj <- coxstuff(x, y, ic, offset = offset)
        }
        nf <- coxstuff.obj$nf
        fail.times <- coxstuff.obj$fail.times
        s <- coxstuff.obj$s
        d <- coxstuff.obj$d
        dd <- coxstuff.obj$dd
        nn <- coxstuff.obj$nn
        nno <- coxstuff.obj$nno

x2<- x^2
oo <- (1.:n)[y >= fail.times[1] ]
sx<-(1/nno[1])*rowSums(x[, oo] * exp(offset[oo]))
s<-(1/nno[1])*rowSums(x2[, oo] * exp(offset[oo]))
w <-  d[1] * (s - sx * sx)


       for(i in 2.:nf) {
           oo <- (1.:n)[y >= fail.times[i-1] & y < fail.times[i] ]
      sx<-(1/nno[i])*(nno[i-1]*sx-rowSums(x[, oo,drop=F] * exp(offset[oo])))
         s<-(1/nno[i])*(nno[i-1]*s-rowSums(x2[, oo,drop=F] * exp(offset[oo])))
       w <- w + d[i] * (s - sx * sx)
        }
        return(w)
}




coxstuff<-
function(x, y, ic, offset = rep(0., length(y)))
{
        fail.times <- unique(y[ic == 1.])
        nf <- length(fail.times)
        n <- length(y)
        nn <- rep(0., nf)
        nno <- rep(0., nf)
        for(i in 1.:nf) {
                nn[i] <- sum(y >= fail.times[i])
                nno[i] <- sum(exp(offset)[y >= fail.times[i]])
        }
        s <- matrix(0., ncol = nf, nrow = nrow(x))
        d <- rep(0., nf)
        #expand d out to a vector of length n
        for(i in 1.:nf) {
                o <- (1.:n)[(y == fail.times[i]) & (ic == 1.)]
                d[i] <- length(o)
}
         oo <- match(y, fail.times)
         oo[ic==0]<-NA
         s<-t(rowsum(t(x),oo))
         s<-s[,-ncol(s)]
        dd <- rep(0., n)
        for(j in 1.:nf) {
                dd[(y == fail.times[j]) & (ic == 1.)] <- d[j]
        }
        return(list(fail.times=fail.times, s=s, d=d, dd=dd, nf=nf, nn=nn, nno=nno))
}


 ocoxvar <-
function(x, y, ic, offset = rep(0., length(y)), coxstuff.obj = NULL)
{
        # computes information elements (var) for cox
        # x is nx by n matrix of expression  values
        nx <- nrow(x)
        n <- length(y)
        yy <- y + (ic == 0.) * (1e-06)
        otag <- order(yy)
        y <- y[otag]
        ic <- ic[otag]
        x <- x[, otag, drop = F]
        offset <- offset[otag]
        if(is.null(coxstuff.obj)) {
                coxstuff.obj <- coxstuff(x, y, ic, offset = offset)
        }
        nf <- coxstuff.obj$nf
        fail.times <- coxstuff.obj$fail.times
        s <- coxstuff.obj$s
        d <- coxstuff.obj$d
        dd <- coxstuff.obj$dd
        nn <- coxstuff.obj$nn
        nno <- coxstuff.obj$nno
        w <- rep(0., nx)
 x2<- x^2
        for(i in 1.:nf) {
                oo <- (1.:n)[y >= fail.times[i]]
      sx<-(1/nno[i])*rowSums(x[, oo] * exp(offset[oo]))
         s<-(1/nno[i])*rowSums(x2[, oo] * exp(offset[oo]))
       w <- w + d[i] * (s - sx * sx)
        }
        return(w)
}


mysvd<-function(x,  n.components=NULL){
# finds PCs of matrix x
  p<-nrow(x)
  n<-ncol(x)

# center the observations (rows)

 x<-t(scale(t(x),center=T,scale=F))

  if(is.null(n.components)){n.components=min(n,p)}
  if(p>n){
    a<-eigen(t(x)%*%x)
    v<-a$vec[,1:n.components,drop=FALSE]
    d<-sqrt(a$val[1: n.components,drop=FALSE])
    
      u<-scale(x%*%v,center=FALSE,scale=d)
 
    
    return(list(u=u,d=d,v=v))
  }
  else{

      junk<-svd(x,LINPACK=TRUE)
      nc=min(ncol(junk$u), n.components)
      return(list(u=junk$u[,1:nc],d=junk$d[1:nc],
                  v=junk$v[,1:nc]))
}
}
superpc.plot.lrtest<- function(object.lrtestcurv){
  plot(object.lrtestcurv$threshold, object.lrtestcurv$lrtest,xlab="Threshold",ylab="Likelihood ratio statistic",type="b")


}
  
superpc.plotcv<-
  
  function (object,  smooth=TRUE, smooth.df=10, ...) 
{

  scor=object$scor.preval

  if(smooth){
    for(j in 1:nrow(scor)){
      scor[j,]=smooth.spline(object$th, scor[j,],df=smooth.df)$y
    }}
  
  ymax=max(scor,qchisq(.95,nrow(scor)))

  matplot(object$th, t(scor), xlab = "Threshold", 
          ylab = "Likelihood ratio test statistic", ylim=c(0,ymax))
  matlines(object$th, t(scor), ...)

  for(j in 1:nrow(scor)){
    abline(h=qchisq(.95,j),lty=2,col=j)
  }

}

superpc.plotshrink.lrtest<- function(object.lrtestred){
  
  n.components<- object.lrtestred$n.components
  
  if(n.components==1) {par(mfrow=c(1,1))}
  
  if(n.components>1) {par(mfrow=c(2,2))}
  
  plot(object.lrtestred$shrinkages, object.lrtestred$lrtest,xlab="Shrinkage amount",ylab="Likelihood ratio test statistic",type="b")
  if(n.components==1){  axis(3,at=object.lrtestred$shrinkages, labels=as.character(object.lrtestred$num.features))}
  
  
  if(n.components>1)
    for(ii in 1:n.components){
      plot(object.lrtestred$shrinkages, object.lrtestred$num.features[,ii],xlab="Shrinkage amount",ylab="Number of features",type="b", log="y")
      title(paste("Component ",as.character(ii),sep=""))
    }
  
}
superpc.predict <- function (object, data, newdata, threshold, n.components=1,  prediction.type=c("continuous", "discrete",
                                                                                  "nonzero"), n.class=2)
{
  
#thresholds the feature scores at "threshold", computes svd based on "data", and then predicts based
#  "newdata"
  
  this.call <- match.call()

  prediction.type <- match.arg(prediction.type)

  which.features <- (abs(object$feature.scores) >= threshold)
  x.sml <- data$x[which.features, ]
  n.pc <- n.components
                                

  x.sml.svd <- mysvd(x.sml, n.components=n.components)



  
  if (prediction.type=="nonzero") {
    if (!is.null(data$featurenames)) {
      out<- data$featurenames[which.features]
    }
    else {
      out<- (1:nrow(data$x))[which.features]
    }
  }

  if (prediction.type=="continuous" | prediction.type=="discrete") {
    cur.v <-scale( t(newdata$x[which.features,]) %*%x.sml.svd$u,  center=FALSE,scale=x.sml.svd$d)

     cur.v0 <-scale( t(data$x[which.features,]) %*%x.sml.svd$u,  center=FALSE,scale=x.sml.svd$d)

#here we obtain the regression coefs of y on the latent factors
# and flip the sign of the factors if the coef is negative

result<-superpc.fit.to.outcome(object, data, cur.v0, print=FALSE)

if(object$type=="survival"){coef=result$coef}
if(object$type=="regression"){coef=result$coef[-1]}



    if (prediction.type=="continuous") {
      out<-scale(cur.v, center=FALSE,scale= sign(coef))
    }
    else if (prediction.type=="discrete") {

     out<-scale(cur.v, center=FALSE,scale= sign(coef))


      for(j in 1:ncol(out)){
        out[,j]<-cut(out[,j],n.class,labels=FALSE)
      }}
  }
  junk <- list(v.pred=out, u = x.sml.svd$u, d = x.sml.svd$d, 
               which.features = which.features, 
               n.components = n.pc, 
               call = this.call, prediction.type=prediction.type)

  return(junk)
}
superpc.predict.red <- function(fit, data, data.test, threshold, n.components=1, n.shrinkage=20, compute.lrtest=TRUE, sign.wt="both"){

  # try reduced predictor on test set

  
  soft.thresh<- function(x,tt){ sign(x)*(abs(x)-tt)*(abs(x)>tt)}


  this.call<- match.call()
  

  type=fit$type

  lrtest.shrink<- rep(NA,n.shrinkage)
  
  cur.vall<- array(NA,c(n.shrinkage,ncol(data$x),n.components))
  cur.vall.test<- array(NA, c(n.shrinkage,ncol(data.test$x),n.components))
  corr.with.full<-matrix(NA,nrow=n.shrinkage, ncol=n.components)
  which.features <- abs(fit$feature.scores) > threshold
  x.sml <- data$x[which.features, ]
  x.svd <- mysvd(x.sml, n.components=n.components)
  cur.v <- scale(t(data$x[which.features, ]) %*%x.svd$u, center=FALSE,scale=x.svd$d)

# flip the sign of the latent factors, if a coef is neg

result<-superpc.fit.to.outcome(fit, data, cur.v, print=FALSE)
if(fit$type=="survival"){coef=result$coef}
if(fit$type=="regression"){coef=result$coef[-1]}

 cur.v<-scale(cur.v, center=FALSE,scale= sign(coef))
##

  sc<-cor(t(data$x), cur.v)

  # don't shrink all of the way to zero 

  maxshrink=max(abs(sc))

  if(sign.wt=="positive"){ maxshrink=max(abs(sc[sc>0]))}
  
  if(sign.wt=="negative"){ maxshrink=max(abs(sc[sc<0]))}

  shrinkages<- seq(0,maxshrink,length=n.shrinkage+1)
  shrinkages= shrinkages[-(n.shrinkage+1)]

  num.features<-matrix(NA,nrow=n.shrinkage, ncol=n.components)
  
  feature.list<-vector("list", n.shrinkage)


  for(i in 1:n.shrinkage){
    cat(i)
    sc2<- soft.thresh(sc,shrinkages[i])
    if(sign.wt=="positive"){sc2[sc2<0]<-0}
    if(sign.wt=="negative"){sc2[sc2>0]<-0}
    nonzero<-sc2!=0
    

    num.features[i,]<- apply(nonzero,2,sum)
    
    junk=vector("list",n.components)
    for(ii in 1:n.components){
      junk[[ii]]<- (1:nrow(data$x))[nonzero[,ii]]
    }
    feature.list[[i]]=junk

    for(ii in 1:n.components){
      cur.vall[i,,ii]<-apply(t(scale(t(data$x[nonzero[,ii],,drop=FALSE]), center=FALSE,scale=1/sc2[nonzero[,ii],ii])),2,sum)
      cur.vall.test[i,,ii]<-apply(t(scale(t(data.test$x[nonzero[,ii],,drop=FALSE]), center=FALSE,scale=1/sc2[nonzero[,ii],ii])),2,sum)
    }}
  cat("",fill=TRUE)

  if(compute.lrtest){ 
    for(i in 1:n.shrinkage){
      if(type=="survival"){
        require(survival)
        junk<- coxph(Surv(data.test$y, data.test$status) ~cur.vall.test[i,,])$loglik
        lrtest.shrink[i]=2*(junk[2]-junk[1])
      }
      else{
        junk<- summary(lm(data.test$y~cur.vall.test[i,,]))
        if(!is.null(junk$fstat)){lrtest.shrink[i]<-junk$fstat[1]}
      }
    }
  }
  
  for(ii in 1:n.components){
    corr.with.full[,ii]=cor(t(cur.vall[,,ii]),cur.v[,ii])
  }

  return(list(shrinkages=shrinkages, lrtest.shrink=lrtest.shrink, corr.with.full=corr.with.full,
              num.features=num.features, feature.list=feature.list, import=sc, v.test=cur.vall.test,
              n.components=n.components, sign.wt=sign.wt, type=type,call=this.call))
  
}



superpc.predict.red.cv <- function(fitred, fitcv, data, threshold, n.shrinkage=30, sign.wt="both"){

 # try reduced predictor on cv folds, via prevalidation

                           
  this.call=match.call()

  type=fitred$type

  n.components=fitred$n.components


  n.fold<-length(fitcv$folds)

  shrinkages<- fitred$shrinkages
  n.shrinkage<-length(shrinkages)
  cur.vall<- array(NA,c(n.shrinkage,ncol(data$x),n.components))

  for(j in 1:n.fold){
    cat(j,fill=TRUE)
    fit.temp<-list(feature.scores=fitcv$featurescores.fold[,j], type=type)
    ii<-fitcv$folds[[j]]
    
    data1<-list(x=data$x[,-ii],y=data$y[-ii],status=data$status[-ii])
    data2<-list(x=data$x[,ii],y=data$y[ii],status=data$status[ii])
    junk<- superpc.predict.red(fit.temp, data1,data2, threshold, n.shrinkage=n.shrinkage, n.components=n.components,compute.lrtest=FALSE, sign.wt=sign.wt)
    cur.vall[,ii,]<-junk$v.test
  }

  lrtest.shrink<-rep(NA,n.shrinkage)

  for(i in 1:n.shrinkage){
    if(type=="survival"){
      require(survival)
      junk<- coxph(Surv(data$y, data$status) ~cur.vall[i,,])$loglik
      lrtest.shrink[i]=2*(junk[2]-junk[1])
    }
    else{
      junk<- summary(lm(data$y~cur.vall[i,,]))
      if(!is.null(junk$fstat)){lrtest.shrink[i]<-junk$fstat[1]}
    }

  }


  return(list(shrinkages=shrinkages, lrtest.shrink=lrtest.shrink, num.features=fitred$num.features, n.components=n.components, v.preval.red=cur.vall, sign.wt=sign.wt, type=type,call=this.call))
}
superpc.train<-
  function (data, type=c("survival","regression")){
    
# computes feature scores for supervised pc analysis
    
  this.call <- match.call()
 type <- match.arg(type)

  
  
  if (is.null(data$status) & type=="survival") {
 stop("Error: survival specified but censoring status is null")
  }
  
  if (type=="survival") {
    feature.scores <- coxfunc(data$x, data$y, data$status)$tt
   }
  else {
    feature.scores <- cor.func(data$x, data$y)$tt
  }

  
  junk <- list( feature.scores=feature.scores, 
               type=type, 
               call = this.call)

  
  class(junk) <- "superpc"
  return(junk)

}

