R/sparseDecom2boot.R

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

#' Convenience wrapper for 2-view eigenanatomy decomposition w/bootstrap
#' initialization.
#'
#' Decomposes two matrices into paired sparse eigenevectors to maximize
#' canonical correlation.
#'
#'
#' @param inmatrix input as inmatrix=list(mat1,mat2). n by p input matrix and n
#' by q input matrix , spatial variable lies along columns.
#' @param inmask optional pair of antsImage masks
#' @param sparseness a c(.,.) pair of values e.g c(0.01,0.1) enforces an
#' unsigned 99 percent and 90 percent sparse solution for each respective view
#' @param nvecs number of eigenvector pairs
#' @param its number of iterations, 10 or 20 usually sufficient
#' @param cthresh cluster threshold pair
#' @param statdir temporary directory if you want to look at full output
#' @param perms number of permutations
#' @param uselong enforce solutions of both views to be the same - requires
#'  matrices to be the same size
#' @param z subject space (low-dimensional space) sparseness value
#' @param smooth smooth the data (only available when mask is used)
#' @param robust rank transform input matrices
#' @param mycoption enforce 1 - spatial orthogonality, 2 - low-dimensional
#' orthogonality or 0 - both
#' @param initializationList initialization for first view
#' @param initializationList2 initialization for 2nd view
#' @param ell1 gradient descent parameter, if negative then l0 otherwise use l1
#' @param nboot n bootstrap runs
#' @param nsamp number of samples e.g. 0.9 indicates 90 percent of data
#' @param doseg boolean to control matrix orthogonality during bootstrap
#' @param priorWeight Scalar value weight on prior between 0 (prior is weak)
#' and 1 (prior is strong).  Only engaged if initialization is used
#' @param verbose activates verbose output to screen
#' @param estimateSparseness effect size to estimate sparseness per vector
#' @return outputs a decomposition of a pair of matrices
#' @author Avants BB
#' @examples
#'
#' \dontrun{
#' mat<-replicate(100, rnorm(20))
#' mat2<-replicate(100, rnorm(20))
#' mydecom<-sparseDecom2boot( inmatrix=list(mat,mat2),
#'   sparseness=c(0.1,0.3) , nvecs=3, its=3, perms=0)
#' wt<-0.666
#' mat3<-mat*wt+mat2*(1-wt)
#' mydecom<-sparseDecom2boot( inmatrix=list(mat,mat3),
#'   sparseness=c(0.2,0.2), nvecs=5, its=10, perms=200 )
#'
#' }
#'
#' @export sparseDecom2boot
sparseDecom2boot <- function(inmatrix, inmask = c(NULL, NULL),
  sparseness = c(0.01, 0.01),
  nvecs = 50, its = 5, cthresh = c(0, 0), statdir = NA, perms = 0,
  uselong = 0,
  z = 0, smooth = 0, robust = 0, mycoption = 1,
  initializationList = list(), initializationList2 = list(),
  ell1 = 0.05, nboot = 10, nsamp = 1, doseg = FALSE,
  priorWeight = 0.0, verbose=FALSE,
  estimateSparseness = 0.2 ) {
  numargs <- nargs()
  if (numargs < 1 | missing(inmatrix)) {
    print(args(sparseDecom2boot))
    return(0)
  }
  nsubj <- nrow(inmatrix[[1]])
  mysize <- round(nsamp * nsubj)
  mat1 <- inmatrix[[1]]
  mat2 <- inmatrix[[2]]
  mymask <- inmask
  cca1out <- 0
  cca2out <- 0
  cca1outAuto <- 0
  cca2outAuto <- 0
  bootccalist1 <- list()
  bootccalist2 <- list()
  nsubs = nrow(mat1)
  allmat1 = matrix( ncol=nvecs*nboot, nrow=ncol(mat1) )
  allmat2 = matrix( ncol=nvecs*nboot, nrow=ncol(mat2) )
  for (i in 1:nvecs) {
    makemat <- matrix(rep(0, nboot * ncol(mat1)), ncol = ncol(mat1))
    bootccalist1 <- lappend(bootccalist1, makemat)
    makemat <- matrix(rep(0, nboot * ncol(mat2)), ncol = ncol(mat2))
    bootccalist2 <- lappend(bootccalist2, makemat)
  }
  if (nsamp >= 0.999999999)
    doreplace <- TRUE else doreplace <- FALSE
  for (boots in 1:nboot) {
    mysample <- sample(1:nsubj, size = mysize, replace = doreplace)
    submat1 <- mat1[mysample, ]
    submat2 <- mat2[mysample, ]
    sublist <- list(submat1, submat2)
#    print(paste("boot", boots, "sample", mysize))
    (myres <- sparseDecom2(
      inmatrix = sublist,
      inmask = mymask,
      sparseness = sparseness,
      nvecs = nvecs,
      its = its,
      cthresh = cthresh,
      perms = 0,
      uselong = uselong,
      z = z,
      smooth = smooth,
      robust = robust,
      mycoption = mycoption,
      initializationList = initializationList,
      initializationList2 = initializationList2,
      ell1 = ell1,
      verbose = verbose ))
    myressum <- abs(diag(cor(myres$projections, myres$projections2)))
    cca1 <- (myres$eig1)
    cca2 <- (myres$eig2)
    if (boots > 1 & TRUE ) {  # try to sort results
      cca1copy <- cca1
      mymult <- matrix(rep(0, ncol(cca1) * ncol(cca1)), ncol = ncol(cca1))
      for (j in 1:ncol(cca1out)) {
        for (k in 1:ncol(cca1)) {
          temp1 <- abs(cca1out[, j])
          temp2 <- abs(cca1[, k])
          mymult[j, k] <- .cosineDist(temp1, temp2)
          # sum( abs( temp1/sum(temp1) - temp2/sum(temp2) ) ) mymult[j,k]<-( -1.0 * cor(
          # temp1, temp2 ) )
        }
      }
      for (ct in 1:(ncol(cca1))) {
        arrind <- which(mymult == min(mymult), arr.ind = T)
        cca1copy[, arrind[1]] <- cca1[, arrind[2]]
        mymult[arrind[1], ] <- 0
        mymult[, arrind[2]] <- 0
      }
      cca1 <- cca1copy
      ###### nextview######
      cca2copy <- cca2
      mymult <- matrix(rep(0, ncol(cca2) * ncol(cca2)), ncol = ncol(cca2))
      for (j in 1:ncol(cca2out)) {
        for (k in 1:ncol(cca2)) {
          temp1 <- abs(cca2out[, j])
          temp2 <- abs(cca2[, k])
          mymult[j, k] <- .cosineDist(temp1, temp2)
          # mymult[j,k]<-sum( abs( temp1/sum(temp1) - temp2/sum(temp2) ) ) mymult[j,k]<-(
          # -1.0 * cor( temp1, temp2 ) )
        }
      }
      for (ct in 1:(ncol(cca2))) {
        arrind <- which(mymult == min(mymult), arr.ind = T)
        cca2copy[, arrind[1]] <- cca2[, arrind[2]]
        mymult[arrind[1], ] <- 0
        mymult[, arrind[2]] <- 0
      }
      cca2 <- cca2copy
    } # end try to sort results
    cca1out <- cca1out + (cca1) # * myressum
    cca2out <- cca2out + (cca2) # * myressum
    bootInds = (( boots - 1 )*nvecs+1):(boots*nvecs)
    allmat1[  ,  bootInds ] = (cca1)
    allmat2[  ,  bootInds ] = (cca2)
    for (nv in 1:nvecs) {
      # if (sparseness[1] > 0)
        bootccalist1[[nv]][boots, ] <- (cca1[, nv])
        # else bootccalist1[[nv]][boots, ] <- (cca1[, nv])
      # if (sparseness[2] > 0)
        bootccalist2[[nv]][boots, ] <- (cca2[, nv])
        # else bootccalist2[[nv]][boots, ] <- (cca2[, nv])
    }
  }

  if ( doseg )
  for ( k in 1:nvecs )
    {
    vec1 = cca1out[,k]
    svec1 = sign( vec1 )
    vec1pos = vec1
    vec1neg = vec1
    vec1pos[ svec1 < 0 ] = 0
    vec1neg[ svec1 > 0 ] = 0
    if ( sum( vec1pos * vec1 ) < sum( vec1neg * vec1 ) ) vec1 = vec1 * (-1)
    if ( sparseness[1] < 0 )
      {
      temp = .eanatsparsify( abs(vec1), abs(sparseness[1]) )
      cca1out[,k] = vec1 * sign( temp )
      } else cca1out[,k] = .eanatsparsify( vec1, sparseness[1] )

      vec2 = cca2out[,k]
      svec2 = sign( vec2 )
      vec2pos = vec2
      vec2neg = vec2
      vec2pos[ svec2 < 0 ] = 0
      vec2neg[ svec2 > 0 ] = 0
      if ( sum( vec2pos * vec2 ) < sum( vec2neg * vec2 ) ) vec2 = vec2 * (-1)
      if ( sparseness[2] < 0 )
        {
        temp = .eanatsparsify( abs(vec2), abs(sparseness[2]) )
        cca2out[,k] = vec2 * sign( temp )
        } else cca2out[,k] = .eanatsparsify( vec2, sparseness[2] )
    }
  if ( estimateSparseness > 1.e-9 )
    {
    sparseness[1]=0
    sparseness[2]=0
    for ( k in 1:nvecs )
      {
      vec1 = colMeans( bootccalist1[[k]] )
      vec2 = colMeans( bootccalist2[[k]] )
      v1sd = apply( bootccalist1[[k]] , FUN=sd, MARGIN=2 )
      v1sd[ v1sd == 0 ] = 1
      v2sd = apply( bootccalist2[[k]] , FUN=sd, MARGIN=2 )
      v2sd[ v2sd == 0 ] = 1
      vec1 = vec1 / v1sd
      vec2 = vec2 / v2sd
      vec1[ abs(vec1) < estimateSparseness ] = 0
      vec2[ abs(vec2) < estimateSparseness ] = 0
      nz1 = sum(abs(sign(vec1)))/length(vec1)
      nz2 = sum(abs(sign(vec2)))/length(vec2)
      if ( verbose ) print( paste( "EstimatedSpar", nz1, nz2) )
      if ( nz1 == 0 ) vec1 = rnorm( length( vec1 ) )
      if ( nz2 == 0 ) vec2 = rnorm( length( vec2 ) )
      cca1out[,k] = vec1
      cca2out[,k] = vec2
      }
    }
  init1 = initializeEigenanatomy( t( cca1out ), inmask[[1]] )
  init2 = initializeEigenanatomy( t( cca2out ), inmask[[2]] )
  ccaout = sparseDecom2(
    inmatrix = inmatrix,
    inmask=c( init1$mask, init2$mask ),
    sparseness = sparseness,
    nvecs = nvecs,
    its = its,
    cthresh = cthresh,
    perms = perms,
    uselong = uselong,
    z = z,
    smooth = smooth,
    robust = robust,
    mycoption = mycoption,
    initializationList=init1$initlist,
    initializationList2=init2$initlist,
    ell1 = ell1,
    priorWeight=priorWeight,
    verbose=verbose )
  jh=matrix( 0, nrow=ncol(inmatrix[[1]]), ncol=ncol(inmatrix[[2]]) )
  colnames(jh)=colnames(inmatrix[[2]])
  rownames(jh)=colnames(inmatrix[[1]])
  for ( i in 1:ncol(allmat1) )
    {
    wh1=which( abs( allmat1[,i] ) > 1.e-10 )
    wh2=which( abs( allmat2[,i] ) > 1.e-10 )
    jh[ wh1, wh2 ] = jh[ wh1, wh2 ] + 1
    }
  return(
    list(
      bootsccan=ccaout,
      eig1=cca1out,
      eig2=cca2out,
      bootccalist1=bootccalist1,
      bootccalist2=bootccalist2,
      allmat1=allmat1,
      allmat2=allmat2,
      init1=init1,
      init2=init2,
      jh=jh )
    )
##### old implementation below #####
  cca1outAuto <- cca1out
  cca2outAuto <- cca2out
  for (nv in 1:nvecs) {
    bootmat <- bootccalist1[[nv]]
    vec1 <- apply(bootmat, FUN = mean, MARGIN = 2)
    # mysd1 <- apply(bootmat,FUN=sd,MARGIN=2) vec1[ mysd1 > 0 ] <- vec1[ mysd1 > 0 ]
    # / mysd1[ mysd1 > 0 ] vec1 <- apply(bootmat,FUN=t.test,MARGIN=2) vec1 <-
    # as.numeric( do.call(rbind, vec1)[,1] )
    vec1[is.na(vec1)] <- 0
    if ((nv > 1) & (abs(sparseness[1] * nvecs) < 1) & TRUE) {
      for (j in 1:(nv - 1)) {
        prevec <- cca1out[, j]
        vec1[prevec > 0] <- 0
      }
      cca1out[, nv] <-.eanatsparsify(vec1, abs(sparseness[1]))
      cca1outAuto[, nv] <- cca1out[, nv] # vec1
    } else {
      cca1out[, nv] <-.eanatsparsify(vec1, abs(sparseness[1]))
      cca1outAuto[, nv] <- cca1out[, nv] # vec1
    }
    ### now vec 2 ###
    bootmat <- bootccalist2[[nv]]
    vec2 <- apply(bootmat, FUN = mean, MARGIN = 2)
    # mysd2 <- apply(bootmat,FUN=sd,MARGIN=2) vec2[ mysd2 > 0 ] <- vec2[ mysd2 > 0 ]
    # / mysd2[ mysd2 > 0 ] vec2 <- apply(bootmat,FUN=t.test,MARGIN=2) vec2 <-
    # as.numeric( do.call(rbind, vec2)[,1] )
    vec2[is.na(vec2)] <- 0
    if ((nv > 1) & (abs(sparseness[2] * nvecs) < 1) & TRUE) {
      for (j in 1:(nv - 1)) {
        prevec <- cca2out[, j]
        vec2[prevec > 0] <- 0
      }
      cca2out[, nv] <-.eanatsparsify(vec2, abs(sparseness[2]))
      cca2outAuto[, nv] <- vec2
    } else {
      cca2out[, nv] <-.eanatsparsify(vec2, abs(sparseness[2]))
      cca2outAuto[, nv] <- vec2
    }
  }
  for (i in 1:ncol(cca1outAuto)) {
    mynorm <- sqrt(sum(cca1outAuto[, i] * cca1outAuto[, i]))
    if (mynorm > 0)
      {
      cca1outAuto[, i] <- cca1outAuto[, i]/mynorm
        .eanatsparsify( cca1outAuto[, i]/mynorm, abs(sparseness[1]) )
      print(sum())
      }
    mynorm <- sqrt(sum(cca2outAuto[, i] * cca2outAuto[, i]))
    if (mynorm > 0)
      {
      cca2outAuto[, i] <-
        .eanatsparsify( cca2outAuto[, i]/mynorm, abs(sparseness[2]))
      }
  }
  fakemask1 <- makeImage(c(1, 1, ncol(mat1)), 1)
  fakemask2 <- makeImage(c(1, 1, ncol(mat2)), 1)
  usefakemask <- c((length(dim(myres$eig1)) == 2), (length(dim(myres$eig2)) ==
    2))
  locmask <- inmask
  if (doseg)
    cca1outAuto <- .matrixSeg(t(cca1outAuto))
  if (dim(cca1outAuto)[2] != nvecs)
    cca1outAuto <- t(cca1outAuto)
  cca1out <- (cca1outAuto)
  if (doseg)
    cca2outAuto <- .matrixSeg(t(cca2outAuto))
  if (dim(cca2outAuto)[2] != nvecs)
    cca2outAuto <- t(cca2outAuto)
  cca2out <- (cca2outAuto)
  ####################################################################################
  if (usefakemask[1])
    init1 <- matrixToImages(t(cca1out), fakemask1) else init1 <- matrixToImages(t(cca1out), locmask[[1]])
  if (usefakemask[2])
    init2 <- matrixToImages(t(cca2out), fakemask2) else init2 <- matrixToImages(t(cca2out), locmask[[2]])
  print("Get Final Results")
  if (!usefakemask[1] & !usefakemask[2])
    maskinit <- locmask
  if (usefakemask[1] & usefakemask[2])
    maskinit <- c(fakemask1, fakemask2)
  if (usefakemask[1] & !usefakemask[2])
    maskinit <- c(fakemask1, locmask[[2]])
  if (!usefakemask[1] & usefakemask[2])
    maskinit <- c(locmask[[1]], fakemask2)
  myres <- sparseDecom2(
    inmatrix = inmatrix,
    inmask = maskinit,
    sparseness = sparseness,
    nvecs = nvecs,
    its = its,
    cthresh = cthresh,
    perms = perms,
    uselong = uselong,
    z = z,
    smooth = smooth,
    robust = robust,
    mycoption = mycoption,
    initializationList = init1,
    initializationList2 = init2,
    ell1 = ell1 )
  ########################################################################
  return(
    list(
      projections = myres$projections,
      projections2 = myres$projections2,
      eig1 = myres$eig1,
      eig2 = myres$eig2,
      ccasummary = myres$ccasummary,
      bootccalist1 = bootccalist1,
      bootccalist2 = bootccalist2,
      cca1outAuto = (cca1outAuto),
      cca2outAuto = (cca2outAuto)
      )
    )
}

.cosineDist <- function(xin, yin) {
  x <- t(as.matrix(xin))
  y <- t(as.matrix(yin))
  return(as.numeric(1 - x %*% t(y)/(sqrt(rowSums(x^2) %*% t(rowSums(y^2))))))
}