.packageName <- "plsgenomics"
### TFA.estimate.R  (2005-04-11)
###
###     Prediction of Transcription Factor Activities using PLS
###
### Copyright 2004-04 Anne-Laure Boulesteix and Korbinian Strimmer
###
### 
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



TFA.estimate<-function(CONNECdata,GEdata,ncomp=NULL,nruncv=0,alpha=2/3,unit.weights=TRUE)
{
n<-nrow(GEdata)
m<-ncol(GEdata)
p<-ncol(CONNECdata)
if (nruncv==0&length(ncomp)>1)
 stop("Since length(ncomp)>1, nruncv must be >0")

if (is.null(ncomp))
 {
 ncomp<-min(n,p)
 }
if (n!=nrow(CONNECdata))
 stop("The number of genes must be the same in the gene expression data and in the CONNEC data")
 
X<-scale(CONNECdata,center=TRUE,scale=TRUE)
X[is.na(X)]<-0
Y<-scale(GEdata,center=TRUE,scale=TRUE)
if (nruncv>0)
 {
 ncomp<-pls.regression.cv(X,Y,ncomp=ncomp,nruncv=nruncv) 
 }

Dx<-diag(1/attributes(X)$"scaled:scale")
Dy<-diag(1/attributes(Y)$"scaled:scale")
pls.out<-pls.regression(X,Y,ncomp=ncomp,Xtest=NULL,unit.weights=unit.weights)
TFA<-Dx%*%pls.out$B%*%solve(Dy)
metafactor<-pls.out$Q

return(list(TFA=TFA,metafactor=metafactor,ncomp=ncomp))
}

### pls.lda.R  (2005-04-06)
###
###     Classification with PLS Dimension Reduction and Linear Discriminan Analysis
###
### Copyright 2004-04 Anne-Laure Boulesteix and Korbinian Strimmer
###
### 
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

pls.lda<-function(Xtrain, Ytrain, Xtest=NULL, ncomp, nruncv=0, alpha=2/3, priors=NULL)
{
ntrain<-nrow(Xtrain)
Ytrain<-as.factor(Ytrain)

if (is.vector(Xtest))
 {
 Xtest<-matrix(Xtest,1,length(Xtest))
 }
if (is.null(Xtest))
 {
 Xtest<-Xtrain
 }
if (nruncv==0&length(ncomp)>1) 
 stop("Since length(ncomp)>1, nruncv must be >0")
 
if (nruncv>0)
 {
 ncomp<-pls.lda.cv(Xtrain,Ytrain,ncomp=ncomp,nruncv=nruncv,alpha=alpha,priors=priors)
 }

pls.out<-pls.regression(Xtrain=Xtrain,Ytrain=transformy(Ytrain),Xtest=NULL,ncomp=ncomp)

Ztrain<-as.data.frame(matrix(pls.out$T,ntrain,ncomp))
Ztrain$y<-Ytrain
Ztest<-as.data.frame(scale(Xtest,center=pls.out$meanX,scale=FALSE)%*%pls.out$R)
lda.out<-lda(formula=y~.,data=Ztrain,priors=priors)
predclass<-predict(object=lda.out,newdata=Ztest)$class

return(list(predclass=predclass,ncomp=ncomp))
}

############################

transformy<-function(y)
{
y<-as.numeric(y)
K<-max(y)
if (K>2)
 {
 Y<-matrix(0,length(y),K)
 for (k in 1:K)
  {
  Y[,k]<-as.numeric(y==k)
  Y[,k]<-Y[,k]-mean(Y[,k])
  }
 }
else
 {
 Y<-matrix(y-mean(y),length(y),1)
 }
 
Y
}



### pls.lda.cv.R  (2005-04-06)
###
###     Determination of the number of latent components to be used for Classification with PLS Dimension Reduction and Linear Discriminant Analysis
###
### Copyright 2004-04 Anne-Laure Boulesteix and Korbinian Strimmer
###
### 
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



pls.lda.cv<-function(Xtrain,Ytrain,ncomp,nruncv=20,alpha=2/3,priors=NULL)
{

if (length(ncomp)==1)
 {
 if (ncomp==1)
  return(ncomp)
 else
  ncomp<-1:ncomp
 }
n<-nrow(Xtrain)
ntrain<-floor(n*2/3)
samp<-matrix(0,ntrain,nruncv)
for (i in 1:nruncv)
 {
 samp[,i]<-sample(n,ntrain)
 }
samp<-as.data.frame(samp)
errorcv<-sapply(samp,FUN=pls.lda.sample,Xtrain,Ytrain,ncomp=ncomp,priors=priors)
meanerror<-apply(errorcv,MARGIN=1,FUN=mean)
ncomp<-ncomp[which.min(meanerror)]
return(ncomp)
}


#####################################

pls.lda.sample<-function(samp,X,Y,ncomp,priors=NULL)
{
errorcv<-numeric(length(ncomp))

for (j in ncomp)
  {
  pls.lda.out<-pls.lda(Xtrain=X[samp,],Ytrain=Y[samp],Xtest=X[-samp,],ncomp=j,nruncv=0,priors=priors)
  errorcv[j]<-sum(pls.lda.out$predclass!=Y[-samp])
  }
errorcv
}
### pls.regression.R  (2005-05-10)
###
###     Multivariate Partial Least Squares Regression
###
### Copyright 2004-05 Anne-Laure Boulesteix and Korbinian Strimmer
###
### Part of the code was adopted from the pls.pcr package by Ron Wehrens
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA





#
# original simpls function from "pls.pcr" package
# modified so that weights R are also returned
# 

standard.simpls <- function (Xtrain, Ytrain, Xtest=NULL, ncomp=NULL)
{
    X<-scale(Xtrain,center=TRUE,scale=FALSE)
    meanX<-attributes(X)$"scaled:center"
    if (is.vector(Ytrain))
     {
     Ytrain<-matrix(Ytrain,length(Ytrain),1)
     }
    Y<-scale(Ytrain,center=TRUE,scale=FALSE)
    n <- dim(X)[1]
    p <- dim(X)[2]
    m <- dim(Y)[2]
    if (is.null(ncomp))
     {
     ncomp<-min(n,p)
     }
     
    S <- crossprod(X, Y)
    RR <- matrix(0, ncol = max(ncomp), nrow = p)
    PP <- matrix(0, ncol = max(ncomp), nrow = p)
    QQ <- matrix(0, ncol = max(ncomp), nrow = m)
    TT <- matrix(0, ncol = max(ncomp), nrow = n)
    VV <- matrix(0, ncol = max(ncomp), nrow = p)
    UU <- matrix(0, ncol = max(ncomp), nrow = n)
    B <- array(0, c(dim(X)[2], dim(Y)[2], length(ncomp)))
    if (!is.null(Xtest))
        {
	if (is.vector(Xtest))
	 {
	 Xtest<-matrix(Xtest,1,length(Xtest))
	 }
        Ypred <- array(0, c(dim(Xtest)[1], m, length(ncomp)))
	}
    for (a in 1:max(ncomp)) {
        qq <- svd(S)$v[, 1]
        rr <- S %*% qq
        tt <- scale(X %*% rr, scale = FALSE)
        tnorm <- sqrt(sum(tt * tt))
        tt <- tt/tnorm
        rr <- rr/tnorm
        pp <- crossprod(X, tt)
        qq <- crossprod(Y, tt)
        uu <- Y %*% qq
        vv <- pp
        if (a > 1) {
            vv <- vv - VV %*% crossprod(VV, pp)
            uu <- uu - TT %*% crossprod(TT, uu)
        }
        vv <- vv/sqrt(sum(vv * vv))
        S <- S - vv %*% crossprod(vv, S)
        RR[, a] <- rr
        TT[, a] <- tt
        PP[, a] <- pp
        QQ[, a] <- qq
        VV[, a] <- vv
        UU[, a] <- uu
        if (!is.na(i <- match(a, ncomp))) {
            B[, , i] <- RR[, 1:a, drop = FALSE] %*% t(QQ[, 1:a,
                drop = FALSE])
            if (!is.null(Xtest))
	     {
	     Xtest<-scale(Xtest,scale=FALSE,center=meanX)
             Ypred[, , i] <- Xtest %*% B[, , i]
	     }
        }
    }
    if (length(ncomp)==1)
     {
     B<-B[,,1]
     }
    if (!is.null(Xtest))
        list(B = B, Ypred = Ypred, P = PP, Q = QQ, T = TT, R=RR, meanX=meanX)
    else list(B = B, P = PP, Q = QQ, T = TT, R=RR, meanX=meanX)
}


# 
# original simpls function returns orthonormal Xscores                                      
# but the weight vectors r_i are NOT standardized to length 1.
#
# instead, this functions returns orthogonal Xscores
# with weight vectors r_i of length 1
#                             

unitr.simpls <- function (Xtrain, Ytrain, Xtest=NULL, ncomp=NULL)
{
  pls.out <- standard.simpls(Xtrain=Xtrain, Ytrain=Ytrain, Xtest=Xtest, ncomp)

  # norm of vector
  euclidian.norm <- function(xvec)
  {
    return( sqrt(sum(xvec*xvec)) )
  }

  # Compute norm of weight vectors r_i
  R.norm <- apply(pls.out$R, 2, euclidian.norm)
  
 
  # Scaling matrices
  if (length(R.norm)==1)
   {
   M<-matrix(R.norm,1,1)
   Mi<-matrix(1/R.norm,1,1)
   }
  else
   {  
   M <- diag(R.norm)
   Mi <- diag(1/R.norm)
   }
  # Transform output matrices
  Rnew <- pls.out$R %*% Mi
  Tnew <- pls.out$T %*% Mi
  
  Qnew <- pls.out$Q %*% M
  Pnew <- pls.out$P %*% M

  # B and Ypred are invariant !


  if (!is.null(Xtest))
        list(B = pls.out$B, Ypred = pls.out$Ypred, P = Pnew, Q = Qnew,
             T = Tnew, R=Rnew, meanX=pls.out$meanX)
    else list(B = pls.out$B, P = Pnew, Q = Qnew, T = Tnew, R=Rnew,
    meanX=pls.out$meanX)
}

######################

pls.regression <- function(Xtrain, Ytrain, Xtest=NULL, ncomp=NULL,  unit.weights=TRUE)
{
  if (unit.weights==TRUE)
  {
    return( unitr.simpls(Xtrain, Ytrain, Xtest=Xtest, ncomp) )
  }
  else
  {
    return( standard.simpls(Xtrain, Ytrain, Xtest=Xtest, ncomp) )
  }
}

### pls.regression.cv.R  (2005-04-06)
###
###    Determination of the number of latent components to be used in PLS regression
###
### Copyright 2004-04 Anne-Laure Boulesteix and Korbinian Strimmer
###
### 
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



pls.regression.cv<-function(Xtrain,Ytrain,ncomp,nruncv=20,alpha=2/3)
{
n<-nrow(Xtrain)
ntrain<-floor(n*2/3)
if (is.vector(Ytrain))
 {
 Ytrain<-matrix(Ytrain,ntrain,1)
 }
samp<-matrix(0,ntrain,nruncv)


if (length(ncomp)==1)
 {
 if (ncomp==1)
  return(ncomp)
 else 
 ncomp<-1:ncomp
 }
  
for (i in 1:nruncv)
 {
 samp[,i]<-sample(n,ntrain)
 }
samp<-as.data.frame(samp)

errorcv<-sapply(samp,FUN=pls.regression.sample,Xtrain,Ytrain,ncomp)
meanerror<-apply(errorcv,MARGIN=1,FUN=mean)

ncomp<-ncomp[which.min(meanerror)]

return(ncomp)
}


#####################################

pls.regression.sample<-function(samp,X,Y,ncomp)
{
lncomp<-length(ncomp)
errorcv<-numeric(lncomp)
pls.out<-pls.regression(X[samp,],Y[samp,],ncomp=ncomp,Xtest=X[-samp,])

for (j in 1:lncomp)
  {
  errorcv[j]<-sum((pls.out$Ypred[,,j]-Y[-samp,])^2)
  }

errorcv
}
### variable.selection.R  (2005-04-06)
###
###     Variable selection using the PLS weights
###
### Copyright 2004-04 Anne-Laure Boulesteix and Korbinian Strimmer
###
### Part of the code was adopted from the pls.pcr package by Ron Wehrens
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

variable.selection<-function(X,Y,nvar=NULL)
{
Xscaled<-scale(X,center=TRUE,scale=TRUE)
Y<-as.numeric(Y)
Y<-Y-mean(Y)
a<-t(Xscaled)%*%Y
if (is.null(nvar))
 {
 nvar<-ncol(X)
 }
return(order(-abs(a))[1:nvar])
}
