## This work is licensed under the Creative Commons
## Attribution-ShareAlike 3.0 Unported License. To view a copy of this
## license, visit http://creativecommons.org/licenses/by-sa/3.0/ or send
## a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain
## View, California, 94041, USA.

## Copyright (c) 2011: Thomas Rusch and Norbert Walchofer

spca_approxOLD <- function(Sigma=NULL,dat=NULL,k,scale=FALSE)
  #@Input:
  # Sigma ... covariance matrix (p times p);
  # A... centered data matrix 
  # k ... max cardinality to check
  #@Output:
  # list with components
  # var... vector of variances for target cardinality
  # rhobreaks... corresponding rho penalties
  # index ... indices with the selected components
  # pattern... the sparse pattern
  {
    if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(A)%*%A
    if(!is.null(Sigma)&&is.null(dat)) dat <- chol(Sigma)
    if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.")
    
    #initialize objects
    rhos <- rep(NA,k)
    p <- dim(A)[2]
    var.out <- rep(NA,k)
    eig.val <- rep(NA,k)
    patt <- rep(0,p)

    #preprocessing
    vars <- rowSums(Sigma*diag(p))
    A <- dat
    
    #initialize algorithm
    res <- which(vars==max(vars))
    x <- dat[,res]/sqrt(sum(dat[,res]^2))
    var.out[1] <- max(vars)
    eig.val[1] <- max(vars)
    rhos[1] <- max(vars)
    patt[res] <- 1 

    if(k==1) return(list(index=res,pattern=patt,variances=eig.val,rhos=rhos))
    for (i in 2:k)
      {
        cat("cardinality = ",i,"\n")
        A[,res] <- 0
        var.tmp <- colSums(x*A)^2
        res.tmp <- which(var.tmp==max(var.tmp))
        rhos[i] <- max(var.tmp)
        res <- c(res,res.tmp)
        mat.tmp <- dat[,res]%*%t(dat[,res])
        eig.tmp <- eigen(mat.tmp)
        x <- eig.tmp$vectors[,1]
        eig.val[i] <- eig.tmp$values[1]
      }
    patt[res] <- 1
    out <- list(index=res,pattern=patt,variances=eig.val,rhos=rhos,subSigma=Sigma[res,res])
    return(out)
  }

spca_approx <- function(Sigma=NULL,dat=NULL,k,scale=FALSE)
  #@Input:
  # Sigma ... covariance matrix (p times p);
  # A... centered data matrix 
  # k ... max cardinality to check
  #@Output:
  # list with components
  # var... vector of variances for target cardinality
  # rhobreaks... corresponding rho penalties
  # index ... indices with the selected components
  # pattern... the sparse pattern
  {
    if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(dat)%*%dat
    if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.")
    if(scale && !is.null(dat)) dat <- scale(dat,scale=FALSE)
    p <- dim(Sigma)[2]
    vars <- rowSums(Sigma*diag(p))
  #  Sigma <- Sigma[,order(vars,decreasing=TRUE)]
  #  Sigma <- Sigma[order(vars,decreasing=TRUE),]
    Sigma <- Sigma[order(vars,decreasing=TRUE),order(vars,decreasing=TRUE)]
  #  if(is.null(dat)) dat <- chol(Sigma)
  #  if(!is.null(dat)) dat <- dat[,order(vars,decreasing=TRUE)]
    ifelse(is.null(dat),dat <- chol(Sigma),dat <- dat[,order(vars,decreasing=TRUE)])
    
    #initialize objects
    rhos <- rep(NA,k)
    var.out <- rep(NA,k)
    eig.val <- rep(NA,k)
    patt <- rep(0,p)

    #preprocessing
    ordi <- order(vars,decreasing=TRUE)
    A <- dat
    
    #initialize algorithm
    #vars <- rowSums(Sigma*diag(p))
    res <- 1
    x <- as.numeric(dat[,res]/sqrt(sum(dat[,res]^2)))
    var.out[1] <- max(vars)
    eig.val[1] <- max(vars)
    rhos[1] <- max(vars)
    patt[res] <- 1 

    if(k==1) return(list(res=res,index=ordi[res],pattern=patt,variances=eig.val,rhos=rhos))

    for (i in 2:k)
      {
#        cat("cardinality = ",i,"\n")
        A[,res] <- 0
        var.tmp <- colSums(x*A)^2
        res.tmp <- which(var.tmp==max(var.tmp))
        rhos[i] <- max(var.tmp)
        res <- c(res,res.tmp)
        mat.tmp <- dat[,res]%*%t(dat[,res])
        eig.tmp <- eigen(mat.tmp)
        x <- eig.tmp$vectors[,1]
        eig.val[i] <- eig.tmp$values[1]
      }
    indi <- ordi[res]
    patt[indi] <- 1
    out <- list(res=res,index=indi,pattern=patt,variances=eig.val,rhos=rhos,Sigma=Sigma,A=dat)
    return(out)
  }

spca_greedy <- function(Sigma=NULL,dat=NULL,k)
  #@Input:
  # Sigma ... covariance or correlation matrix (p times p correlations);
  # A... centered data matrix 
  # k ... max cardinality to check
  #@Output:
  # list with components
  # var... vector of variances for target cardinality
  # rhobreaks... corresponding rho penalties
  # index ... indices with the selcted components
  # pattern... the sparse pattern
  {
    if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(A)%*%A
    if(!is.null(Sigma)&&is.null(dat)) dat <- chol(Sigma)
    if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.")
    
    #initialize objects
    rhos <- rep(NA,k)
    p <- dim(A)[2]
    var.out <- rep(NA,k)
    eig.val <- rep(NA,k)
    patt <- rep(0,p)

    #preprocessing
    vars <- rowSums(Sigma*diag(p))
    A <- dat
    
    #initialize algorithm
    res <- which(vars==max(vars))
    x <- dat[,res]/sqrt(sum(dat[,res]^2))
    var.out[1] <- max(vars)
    eig.val[1] <- max(vars)
    rhos[1] <- max(vars)
    patt[res] <- 1 

    if(k==1) return(list(index=res,pattern=patt,variances=eig.val,rhos=rhos))

    for (i in 2:k)
      {
        A[,res] <- 0
        var.tmp <- colSums(x*A)^2#hier was tun
        res.tmp <- which(var.tmp==max(var.tmp))
        rhos[i] <- max(var.tmp)
        res <- c(res,res.tmp)
        mat.tmp <- dat[,res]%*%t(dat[,res])
        eig.tmp <- eigen(mat.tmp)
        x <- eig.tmp$vectors[,1]
        eig.val[i] <- eig.tmp$values[1]
      }
    patt[res] <- 1
    out <- list(index=res,pattern=patt,variances=eig.val,rhos=rhos)
    return(out)
  }

check_optimality <- function(Sigma,subset)
  {
    if(all(subset==c(1,2,5,3,6))) return(list(optimal=TRUE,rho=2.18e04))
    return(list(optimal=FALSE,rho=NA))
  }