R/sparseDecom.R

In neuroconductor-devel/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Convenience wrapper for eigenanatomy decomposition.
#'
#' Decomposes a matrix into sparse eigenevectors to maximize explained
#' variance. Note: we do not scale the matrices internally.
#' We leave scaling choices to the user.
#'
#' @param inmatrix n by p input images , subjects or time points by row ,
#' spatial variable lies along columns
#' @param inmask optional antsImage mask
#' @param sparseness lower values equal more sparse
#' @param nvecs number of vectors
#' @param its number of iterations
#' @param cthresh cluster threshold
#' @param statdir place on disk to save results
#' @param z u penalty, experimental
#' @param smooth smoothness eg 0.5
#' @param initializationList see initializeEigenanatomy
#' @param mycoption 0, 1 or 2 all produce different output 0 is combination
#'   of 1 (spatial orthogonality) and 2 (subject space orthogonality)
#' @param robust rank transform input data - good for data checking
#' @param ell1 the ell1 grad descent param
#' @param getSmall try to get smallest evecs (bool)
#' @param verbose activates verbose output
#' @param powerit alternative power iteration implementation, faster
#' @param priorWeight scalar weight typically in range zero to two
#' @param maxBased boolean that chooses max-based thresholding
#' @return outputs a decomposition of a population or time series matrix
#' @author Avants BB
#' @examples
#'
#' mat<-replicate(100, rnorm(20))
#' mydecom<-sparseDecom( mat )
#' mat<-scale(mat)
#' mydecom2<-sparseDecom( mat )
#' # params that lead to algorithm similar to NMF
#' mydecom3<-sparseDecom( mat, z=1, sparseness=1 )
#'
#' \dontrun{
#' # for prediction
#' if ( usePkg("randomForest") & usePkg("spls")  & usePkg('BGLR') ) {
#' data(lymphoma) # from spls
#' training<-sample( rep(c(TRUE,FALSE),31)  )
#' sp<-0.02 ; myz<-0
#' ldd<-sparseDecom( lymphoma$x[training,], nvecs=5 , sparseness=( sp ),
#'   mycoption=1, z=myz ) # NMF style
#' traindf<-data.frame( lclass=as.factor(lymphoma$y[ training  ]),
#'   eig = lymphoma$x[training,]  %*% as.matrix(ldd$eigenanatomyimages ))
#' testdf<-data.frame(  lclass=as.factor(lymphoma$y[ !training ]),
#'  eig = lymphoma$x[!training,] %*% as.matrix(ldd$eigenanatomyimages ))
#' myrf<-randomForest( lclass ~ . ,   data=traindf )
#' predlymp<-predict(myrf, newdata=testdf)
#' print(paste('N-errors:',sum(abs( testdf$lclass != predlymp ) ),
#'   ' non-zero ',sum(abs( ldd$eigenanatomyimages ) > 0 ) ) )
#' # compare to http://arxiv.org/pdf/0707.0701v2.pdf
#' # now SNPs
#' data(mice)
#' snps<-quantifySNPs( mice.X, shiftit = TRUE )
#' numericalpheno<-as.matrix( mice.pheno[,c(4,5,13,15) ] )
#' nfolds<-6
#' train<-sample( rep( c(1:nfolds), 1800/nfolds ) )
#' train<-( train < 4 )
#' lrmat<-lowrankRowMatrix( as.matrix( snps[train,] ) ,  50 )
#' lrmat=scale(lrmat)
#' snpd<-sparseDecom( lrmat-min(lrmat), nvecs=20 , sparseness=( 0.001), z=-1 )
#' projmat<-as.matrix( snpd$eig )
#' snpse<-as.matrix( snps[train, ]  ) %*% projmat
#' traindf<-data.frame( bmi=numericalpheno[train,3] , snpse=snpse)
#' snpse<-as.matrix( snps[!train, ]  ) %*% projmat
#' testdf <-data.frame( bmi=numericalpheno[!train,3] , snpse=snpse )
#' myrf<-randomForest( bmi ~ . , data=traindf )
#' preddf<-predict(myrf, newdata=testdf )
#' cor.test(preddf, testdf$bmi )
#' plot(preddf, testdf$bmi )
#' } # check for packages
#' # prior-based example
#' set.seed(123)
#' ref<-antsImageRead( getANTsRData("r16"))
#' ref<-iMath(ref,"Normalize")
#' mi<-antsImageRead( getANTsRData("r27"))
#' mi2<-antsImageRead( getANTsRData("r30"))
#' mi3<-antsImageRead( getANTsRData("r62"))
#' mi4<-antsImageRead( getANTsRData("r64"))
#' mi5<-antsImageRead( getANTsRData("r85"))
#' refmask<-getMask(ref)
#' refmask<-iMath(refmask,"ME",2) # just to speed things up
#' ilist<-list(mi,mi2,mi3,mi4,mi5)
#' for ( i in 1:length(ilist) )
#' {
#' ilist[[i]]<-iMath(ilist[[i]],"Normalize")
#' mytx<-antsRegistration(fixed=ref , moving=ilist[[i]] ,
#'   typeofTransform = c("Affine") )
#' mywarpedimage<-antsApplyTransforms(fixed=ref,moving=ilist[[i]],
#'   transformlist=mytx$fwdtransforms)
#' ilist[[i]]=mywarpedimage
#' }
#' mat=imageListToMatrix( ilist , refmask )
#' kmseg=kmeansSegmentation( ref, 3, refmask )
#' initlist=list()
#' for ( k in 1:3 )
#'  initlist[[k]]=
#'    thresholdImage(kmseg$probabilityimages[[k]],0.1,Inf) *
#'    kmseg$probabilityimages[[k]]
#' eanat<-sparseDecom( mat,
#'   inmask=refmask, ell1=0.1,
#'   sparseness=0.0, smooth=0.5, verbose=1,
#'   initializationList=initlist, cthresh=25,
#'   nvecs=3, priorWeight=0.5 )
#' ee=matrixToImages( eanat$eigenanatomyimages, refmask )
#' eseg=eigSeg(  refmask, ee )
#' priormat=imageListToMatrix( initlist, refmask )
#' cor( t(eanat$eigenanatomyimages), t(priormat) )
#' plot( ref, eseg )
#' }
#' @export sparseDecom
sparseDecom <- function(inmatrix = NA, inmask = NULL,
                        sparseness = 0.1,
                        nvecs = 10,
                        its = 5, cthresh = 50,
                        statdir = NA, z = 0, smooth = 0, initializationList = list(),
                        mycoption = 0, robust = 0, ell1 = 1, getSmall = 0, verbose=FALSE,
                        powerit=0, priorWeight=0,
                        maxBased=FALSE ) {
  numargs <- nargs()
  if (numargs < 1 | missing(inmatrix)) {
    cat(" sparseDecom( inmatrix=NA,  inmask=NULL , sparseness=0.01 , nvecs=50 , its=5 , cthresh=250 ) \n")
    return(0)
  }
  idim=3
  if (!is.null(inmask)) {
    inmask = check_ants(inmask)
    idim = inmask@dimension
    if ( sum( inmask >= 0.5 ) != ncol( inmatrix )  ) {
      stop("mask size not equal to matrix column count")
    }
  } else {
    inmask = new("antsImage", "float", idim)
  }
  verbose = as.numeric( verbose )
  time1 <- (Sys.time())
  if ( robust > 0 )
  {
    outval = .Call( "eigenanatomyCpp",
                    robustMatrixTransform(inmatrix),
                    inmask, sparseness, nvecs, its, cthresh, z, smooth,
                    initializationList, mycoption, ell1, verbose, powerit,
                    priorWeight, maxBased,
                    PACKAGE="ANTsR" )
  } else {
    outval = .Call( "eigenanatomyCpp",
                    inmatrix,
                    inmask, sparseness, nvecs, its, cthresh, z, smooth,
                    initializationList, mycoption, ell1, verbose, powerit,
                    priorWeight, maxBased,
                    PACKAGE="ANTsR" )
  }
  time2 <- (Sys.time())
  outval = lappend( outval,  (time2 - time1) )
  names(outval)[length(outval)]='computationtime'
  if ( verbose )
  {
    temp=lm( inmatrix ~  ( inmatrix %*% t(outval$eigenanatomyimages) ) )
    reconmat = predict( temp )
    reconerr = 0
    for ( i in 1:ncol(inmatrix) )
    {
      temp = abs( cor(inmatrix[,i],reconmat[,i]) )
      reconerr = reconerr + temp
    }
    reconerr = reconerr / ncol( inmatrix )
    #  reconerr=mean( abs( inmatrix - reconmat ) )
    outval[length(outval)-1]=reconerr
    #  names(outval)[length(outval)-1]='meanReconstructionError'
    outval[[length(outval)+1]]=reconmat
    names(outval)[length(outval)]='Reconstruction'
  }
  return( outval )
  #  return(list(projections = mydecomp, eigenanatomyimages = fnl, umatrix = fnu,
}