R/spectralNMF.R

Defines functions spectralNMFList getNMFInputMatrix nonNegativePreprocessing initializeNMFModel scaleNMFResult checkForRedundantSources removeRedundantSources includeRedundantSources upsampleNMFResult completeSourceMatrix predictNNLS computeNMFResidu

Documented in checkForRedundantSources completeSourceMatrix computeNMFResidu getNMFInputMatrix includeRedundantSources initializeNMFModel nonNegativePreprocessing predictNNLS removeRedundantSources scaleNMFResult spectralNMFList upsampleNMFResult

# Project: spectralAnalysis_git
# 
# Author: nsauwen <nicolas.sauwen@openanalytics.eu>
#
# Description: Run NMF analyses on spectral data
#
# Changes:
#  * version 0.3.2 
#      + change of extractor (S4)
#      + object replaces by object for consistency 
#
###############################################################################




#' Perform Non-Negative Matrix factorization on spectral data
#' 
#' @param object \code{\link{SpectraInTime-class}}
#' @param rank number of NMF components to be found
#' @param method name of the NMF method to be used. "PGNMF" (default), "HALSacc" 
#' and "semiNMF" are methods derived from the hNMF package. All methods from the NMF package are also available.
#' @param initSpectralData this can be a list of spectralData objects, containing 
#' the pure component spectra. It can also be either of the NMF factor matrices with initial values 
#' @param nruns number of NMF runs. It is recommended to run the NMF analyses multiple
#' times when random seeding is used, to avoid a suboptimal solution
#' @param subsamplingFactor subsampling factor used during NMF analysis
#' @param checkDivergence Boolean indicating whether divergence checking should be performed
#' @return Scaled NMF model (in accordance with the NMF package definition)
#' @author Nicolas Sauwen
#' @export
spectralNMF        <- function ( object , rank, method = "PGNMF", initSpectralData = NULL, nruns = 10, subsamplingFactor = 3, checkDivergence = TRUE) {
  
  spectra          <- getNMFInputMatrix( object , method)
  seed             <- initializeNMFModel( initSpectralData, spectra = spectra, wavelengths = getWavelengths( object ) )
  # Check if NMF initialization is consistent with specified rank
  if(!is.null(seed)){
	  W0               <- basis(seed)
	  if(ncol(W0) > rank) {
		  stop("Number of provided pure component spectra is larger than specified NMF rank.")    
	  }
	  else if(ncol(W0) < rank){
		  # In that case we complete the initial matrix with random sources
		  warning("Number of specified source components is lower than specified rank. Source matrix will be completed with random vectors.")
		  nruns <- 10
#		  checkDivergence <- FALSE
	  }
  }
  timePointsList <- list()
  timePointsList[[1]] <- getTimePoints( object )
  NMFResult         <- runNMF(spectra, rank, method, seed, nruns, checkDivergence, timePointsList, subsamplingFactor)
  NMFResult         <- NMFResult[[1]]
  NMFResult        <- scaleNMFResult( NMFResult )
  return(NMFResult)
}


#' Perform Non-Negative Matrix factorization on list of SPC files
#' 
#' @param objectList list of SPC files
#' @param rank number of NMF components to be found
#' @param method name of the NMF method to be used, consult the help
#' of the 'nmf' function from the NMF package for the methods available by default
#' @param initSpectralData list of SPC files containing pure component spectra
#' @param nruns number of NMF runs.
#' @param subsamplingFactor subsampling factor used during NMF analysis
#' @param checkDivergence Boolean indicating whether divergence checking should be performed
#' @return list of NMF models
#' @importFrom NMF basis coef nmfModel
#' @author Nicolas Sauwen
#' @export
spectralNMFList             <-  function( objectList ,  rank , method = "PGNMF" , initSpectralData = NULL, nruns = 10, subsamplingFactor = 3, checkDivergence = TRUE) {
  
#  spectra                   <- c() # suggestion better preset object size than growing 
  timePointsList            <- list()
  nObservations 			<- 0
  nWavelengths 				<- ncol(   getSpectra( objectList[[1]] ) )
  
  for( iFile in 1:length(objectList)) {
    nObservations            <-  nObservations + nrow(  getSpectra( objectList[[iFile]] ) )
  }
  
  spectra <- matrix(0, nWavelengths, nObservations)
  count <- 1 # Counter keeps track of where to put the data matrices per experiment in the global spectra matrix
  
  #suggestion seq_along , evem better avoid for loop , use seperate function parts  
  for( iFile in 1:length(objectList)) {
    spectra_temp            <-  getNMFInputMatrix(objectList[[iFile]])
    count_temp 				<- ncol(spectra_temp)
    if(nrow(spectra_temp) != nrow(spectra)) {
      stop("NMF analysis cannot be executed. Please make sure all spectral files cover the same wavelength range.")    
    }
    spectra[ , count:(count+count_temp-1)]	<-  spectra_temp
    timePointsList[[iFile]] 				<- getTimePoints( objectList[[iFile]] )
    count <- count + count_temp
  }
  
  seed             <- initializeNMFModel(initSpectralData, spectra = spectra, wavelengths = getWavelengths( objectList[[1]] ) )
  # Check if NMF initialization is consistent with specified rank
  W0               <- basis(seed)
  if(ncol(W0) > rank) {
    stop("Number of provided pure component spectra is larger than specified NMF rank.")    
  }
  else if(ncol(W0) < rank){
    # In that case we complete the initial matrix with random sources
    warning("Number of specified source components is lower than specified rank. Source matrix will be completed with random vectors.")
    nruns <- 10
    checkDivergence <- FALSE
  }
  NMFResultList        <- runNMF(spectra, rank, method, seed, nruns, checkDivergence, timePointsList, subsamplingFactor)
  
  return(NMFResultList)
}


if( 0 == 1 ) {
  # check NFM function in combination with spectralApps helper functions
  library( spectralApps )
  allSPCData        <-  loadAllobjects(dataDir)
  data              <-  stackSpcData( data , wavelengthSelect = c( 650, 1800 ) , 
      timeSelect = c(1, 5) )
  spectra          <-  data
  rank             <-  2
  method = "PGNMF"
  seed = NULL
  nruns = 10
  checkDivergence = NULL 
  
  # bug spectral NMF 
  
  runNMF( spectra = dataForNMF , rank = NCOMPONENTS , seed = initNmfModel ,
      subsamplingFactor = 3 , timePointsList = list( timePoints  = 1: dim(dataForNMF )[2] )  )
  
  spectra = dataForNMF ; rank = NCOMPONENTS ; seed = initNmfModel ;
  subsamplingFactor = 3 ; timePointsList = list( timePoints  = 1: dim(dataForNMF )[2] ) 
  
  method = "PGNMF";  nruns = 10; checkDivergence = NULL; timePointsList = NULL; subsamplingFactor = 3
  
}

#' Actual NMF analysis
#' 
#' @param spectra spectral input matrix, with wavelengths as its rows and time points as its columns
#' @param rank number of NMF components to be found
#' @param method name of the NMF method to be used, consult the help
#' of the 'nmf' function from the NMF package for the methods available by default
#' @param seed nmfModel object containing initialization of the factor matrices
#' @param nruns number of NMF runs. It is recommended to run the NMF analyses multiple
#' times when random seeding is used, to avoid a suboptimal solution
#' @param checkDivergence Boolean indicating whether divergence checking should be performed, defaults to \code{TRUE}
#' @param timePointsList list of time point vectors of the individual experiments
#' @param subsamplingFactor subsampling factor used during NMF analysis
#' @return Resulting NMF model (in accordance with the NMF package definition)
#' @importFrom NMF basis nmf nmfAlgorithm setNMFMethod
#' @importFrom hNMF HALSacc PGNMF semiNMF
#' @importFrom stats runif
#' @importFrom utils combn
#' @author Nicolas Sauwen
#' @export
runNMF                      <-  function (spectra, rank, method = "PGNMF", seed = NULL, nruns = 10, checkDivergence = TRUE, timePointsList = NULL, subsamplingFactor = 3) {
  
#  # Check if NMF method is already available in the NMF package
#  strComp                  <- match( nmfAlgorithm() , method )
#  if( sum( is.na( strComp ) ) == length(strComp)) {
#    setNMFMethod( method , get( method ) ) # necessary to set method if you provide in in nmf function?
#  }
  nTimePoints <- ncol(spectra)
  
  if(subsamplingFactor != 1) {
    sampleInds <- seq(1,ncol(spectra),subsamplingFactor)
    spectra        <- spectra[, sampleInds]
    if( !is.null( seed ) ) {
      W0           <- basis(seed) 
      H0           <- coef(seed)
      H0           <- H0[, seq(1,ncol(H0),subsamplingFactor), drop = FALSE] 
      seed         <- nmfModel( W = W0 , H = H0 )
    }
  }
  
  if( is.null( seed ) ) {
#    NMFResult               <- nmf( spectra , rank = rank, method = get(method), nrun = nruns, checkDivergence = F, .options="-p")
	
	NMFResult <- NULL
	residu <- Inf
	
	for(iRun in 1:nruns) {	
		W0 <- matrix(runif(rank*nrow(spectra)), nrow = nrow(spectra), ncol = rank)	
		H0 <- matrix(runif(rank*ncol(spectra)), nrow = rank, ncol = ncol(spectra))	
		seed <- nmfModel( W = W0 , H = H0 )
		if(method == "PGNMF"){
			NMFResult_temp <- PGNMF(spectra, nmfMod = seed, checkDivergence = F)
		}
		else if(method == "HALSacc"){
			NMFResult_temp <- HALSacc(spectra, nmfMod = seed, checkDivergence = F)
		}
		else if(method == "semiNMF"){
			NMFResult_temp <- semiNMF(spectra, nmfMod = seed, checkDivergence = F)
		}
		W_temp <- basis(NMFResult_temp)
		H_temp <- coef(NMFResult_temp)
		residu_temp <- norm(spectra - W_temp%*%H_temp,'f')
		if(residu_temp < residu) {
			NMFResult <- NMFResult_temp
			residu <- residu_temp
		}
	}
  }
  else{
#      NMFResult             <-  nmf( spectra , rank = rank , method = get(method) , nrun = 1 , seed = seed , checkDivergence = checkDivergence, .options = "-p" )
	  W0 <- basis(seed)
	  if(ncol(W0) == rank) {
		  if(method == "PGNMF"){
			  NMFResult <- PGNMF(spectra, nmfMod = seed, checkDivergence = checkDivergence)
		  }
		  else if(method == "HALSacc"){
			  NMFResult <- HALSacc(spectra, nmfMod = seed, checkDivergence = checkDivergence)
		  }
		  else if(method == "semiNMF"){
			  NMFResult <- semiNMF(spectra, nmfMod = seed, checkDivergence = checkDivergence)
		  }
	  }
	  else { # This means that the source matrix has to be completed with random vector(s)
		  
		  NMFResult <- NULL
		  residu <- Inf
		  W0_orig <- W0
		  
		  for(iRun in 1:ceiling(nruns/2)) {	
			  W0 <- completeSourceMatrix(W0_orig, rank)
			  seed <- initializeNMFModel( initSpectralData = W0, spectra = spectra)
			  iter <- 1
			  maxIter <- 4
			  overlap <- FALSE # Boolean to indicate when random source vectors start overlapping too much with initialized sources
			  NMFResult_prev <- NULL
			  while(iter <= maxIter & !overlap) { #This repetition is done to be able to keep using checkDivergence = T (although with start initialization with 1 random vector)
				  if(method == "PGNMF"){
					  NMFResult_temp <- PGNMF(spectra, nmfMod = seed, checkDivergence = T)
				  }
				  else if(method == "HALSacc"){
					  NMFResult_temp <- HALSacc(spectra, nmfMod = seed, checkDivergence = T)			  
				  }
				  else if(method == "semiNMF"){
					  NMFResult_temp <- semiNMF(spectra, nmfMod = seed, checkDivergence = T)
				  }
				  NMFResult_temp <- scaleNMFResult(NMFResult_temp)
				  W_temp <- basis(NMFResult_temp)
				  checkOverlapMat <- t(W_temp)%*%W_temp
				  checkOverlapInds <- upper.tri(checkOverlapMat)
				  if(ncol(W0_orig) > 1){
					  skipInds <- t(combn(1:ncol(W0_orig),2))
					  checkOverlapInds[skipInds] <- FALSE
				  }
				  checkOverlap <- checkOverlapMat[checkOverlapInds]
				  # Here we are checking the correlation coefficients between the randomly initialized sources and all sources
				  if(sum(checkOverlap > 0.9) > 0){
					  NMFResult_temp <- NMFResult_prev
					  overlap <- TRUE
				  }
				  W_temp <- cbind(W0_orig, W_temp[,(ncol(W0_orig)+1):ncol(W_temp)])
				  NMFResult_prev <- NMFResult_temp
				  seed <- initializeNMFModel( initSpectralData = W_temp, spectra = spectra)				  
				  iter <- iter + 1
			  }
			  W_temp <- basis(NMFResult_temp)
			  H_temp <- coef(NMFResult_temp)
			  residu_temp <- norm(spectra - W_temp%*%H_temp,'f')
			  if(residu_temp < residu) {
				  NMFResult <- NMFResult_temp
				  residu <- residu_temp
			  }
		  }
	  }
  }
  
  # Split abundances per file and save in a list of NMF model objects 
  W <- basis(NMFResult)
  H <- coef(NMFResult)
  
  if(is.null(timePointsList)) {
    timePointsList <- c(1:nTimePoints)
  }
  
  NMFResultList  <- list()
  ind1           <- 1
  ind2           <- floor((length(timePointsList[[1]])-1)/subsamplingFactor)+1
  currentInd <- 1
  shift1 <- 0
  shift2 <- 0
  
  for(iFile in 1:length(timePointsList)) {
    W_temp       <- W
    H_temp       <- H[, (ind1:ind2), drop = FALSE]
    if(iFile != length(timePointsList)) {
      currentInd <- currentInd + length(timePointsList[[iFile]])
      shiftInd <- sampleInds-currentInd
      shiftInd <- which(shiftInd >= 0)[1]
      shift2 <- sampleInds[shiftInd] - currentInd
      
      ind1       <- ind1 + floor((length(timePointsList[[iFile]])-shift1-1)/subsamplingFactor) + 1
      ind2       <- ind2 + floor((length(timePointsList[[iFile+1]])-shift2-1)/subsamplingFactor) + 1
    }
    NMFResultList[[iFile]] <- nmfModel(W = W_temp, H = H_temp)
    if(subsamplingFactor != 1) {
      NMFResultList[[iFile]]      <- upsampleNMFResult(NMFResultList[[iFile]], timePointsList[[iFile]], subsamplingFactor, shift1)
      NMFResultList[[iFile]]      <- scaleNMFResult(NMFResultList[[iFile]])
    }
    shift1 <- shift2
  }    
  gc()
  return( NMFResultList )
}



#' Extract spectral input matrix from SPC file and condition properly for NMF
#' 
#' @param object object of the 'spectralData' class, such as a raw SPC file
#' @param method name of the NMF method to be used.
#' @return spectral matrix, with wavelengths as its rows and time points as its columns
#' @author Nicolas Sauwen
#' @export
getNMFInputMatrix          <- function(object, method = "") {
  
  spectra                  <-  getSpectra( object )
  spectraPreprocessed      <-  nonNegativePreprocessing( spectra, method )
  return( spectraPreprocessed )
}


# seperate from getNMFInputMatrix because needed in spectralApps

#'condition datamatrix to input in and condition properly for NMF
#' @details put negative values to zero, transpose, an add small value zero row (wavelength with only zeros) 
#' 
#' @param spectra matrix of spectra 
#' @param method name of the NMF method to be used.
#' @return matrix, with wavelengths as its rows and time points as its columns
#' @export
nonNegativePreprocessing    <-  function( spectra, method = "") {
  zeroInds <- which(spectra < 0)
  if(length(zeroInds) > 0 & method != "semiNMF"){
    spectra[ spectra < 0 ]    <-  0
    warning("Input data contain negative values. These were reset to zero for NMF analysis")
  }
  
  spectraT                  <-  t( spectra )
  
  zeroRow                   <-  which( ( apply( spectraT , 1 , sum ) ) == 0 )
  spectraT[ zeroRow , ]     <- 1e-16
  spectraT[ spectraT == 0 ] <- 1e-16
  return(spectraT)
}


if( 0 == 1 ) {
  
# debug not working in app 
  initSpectralData           <-  NULL
  W0Init                     <-  nonNegativePreprocessing( initialization )
  W0 = W0Init  ; spectra = dataForNMF ; wavelengths = dataInput$wavelengths
  initNmfModel               <-  initializeNMFModel( W0 = W0Init  , spectra = dataForNMF ,  wavelengths = dataInput$wavelengths )   
  
}

#' Initialize NMF model with initial spectral data
#' 
#' @param initSpectralData this can be a list of spectralData objects, containing 
#' the pure component spectra. It can also be either of the NMF factor matrices with initial values 
#' @param spectra spectral matrix, with wavelengths as its rows and time points as its columns
#' @param wavelengths vector of wavelength values
#' @importFrom nnls nnls
#' @importFrom stats coefficients approx
#' @importFrom NMF nmfModel
#' @export
initializeNMFModel     <- function(initSpectralData, spectra, wavelengths = NULL) {
  if( !is.list( initSpectralData ) & !is.matrix( initSpectralData ) ){
    NMFInit            <- NULL
    return(NMFInit)
  }  
  else if( is.list( initSpectralData ) ) {
	rank               <-  length(initSpectralData)
    W0                 <-  matrix( 0 , nrow(spectra), rank)
    
    for(iList in 1:rank) {
      wavelengths_0    <-  getWavelengths( initSpectralData[[iList]] )
      spectra_0        <- getSpectra( initSpectralData )[[iList]][1,]
      if(abs(wavelengths[1] - wavelengths_0[1]) > 1e-3 | abs(wavelengths[length(wavelengths)] - wavelengths_0[length(wavelengths_0)]) > 1e-3) {
        stop("Please make sure the pure component SPC files cover the same wavelength range as the input SPC file(s).")    
      }
      if(length(wavelengths) != length(wavelengths_0)) {
        spectra_0      <- approx(wavelengths_0, spectra_0, xout = wavelengths) 
        spectra_0      <- spectra_0$y
      }    
      W0[,iList]       <-  spectra_0
    }
  } 
  else if(is.matrix( initSpectralData ) & nrow(initSpectralData) == nrow(spectra)) {
    W0                 <-  initSpectralData
    rank               <-  ncol( W0 )
  }
  
  else if(is.matrix( initSpectralData ) & ncol(initSpectralData) == ncol(spectra)) {
	  H0                 <-  initSpectralData
	  rank               <-  nrow( H0 )
  }
  else {
	  W0                 <-  initSpectralData
	  rank               <-  ncol( W0 )
  }
  
  if(is.na(match("H0", ls()))){
	  nSignals           <-  ncol(spectra)
	  H0                 <-  matrix( 0 , rank , nSignals)
	  for( iSignal in 1:nSignals ) {
		  nlsFit           <-  nnls::nnls( W0 , spectra[,iSignal] )   
		  H0[,iSignal]     <-  coefficients(nlsFit)
	  }
  }
  
  if(is.na(match("W0", ls()))){
	  nFeatures           <-  nrow(spectra)
	  W0                  <-  matrix( 0 , rank, nFeatures)
	  for( iFeature in 1:nFeatures ) {
		  nlsFit           <-  nnls::nnls( t(H0) , spectra[iFeature, ] )   
		  W0[,iFeature]     <-  coefficients(nlsFit)
	  }
	  W0 <- t(W0)
  }
  
  # Proper scaling of the factor matrices
  N0                 <- sqrt(diag(t(W0)%*%W0))
  N_W0               <- matrix(N0,nrow = nrow(W0),ncol = ncol(W0),byrow = T)
  W0                 <- W0/N_W0
  N_H0               <- matrix(N0,nrow = nrow(H0),ncol = ncol(H0))
  H0                 <- H0*N_H0
  
  NMFInit            <-  nmfModel( W = W0 , H = H0 )
  NMFInit            <-  scaleNMFResult(NMFInit)
  
  return(NMFInit)
}


#' Apply fixed scaling to NMF model matrices by normalizing the basis vectors
#' 
#' @param NMFResult Fitted NMF model
#' @return NMFResult Rescaled NMF model
#' @importFrom NMF basis coef
#' @author Nicolas Sauwen
#' @export
scaleNMFResult      <- function(NMFResult) {
  
  W                 <- basis(NMFResult)
  N                 <- sqrt(diag(t(W)%*%W))
  N_W               <- matrix(N,nrow = nrow(W),ncol = ncol(W),byrow = T)
  W                 <- W/N_W
  H                 <- coef(NMFResult)
  N_H               <- matrix(N,nrow = nrow(H),ncol = ncol(H))
  H                 <- H*N_H
  NMF::basis(NMFResult)  <- W
  NMF::coef(NMFResult)   <- H
  return(NMFResult)
  
}


#' Check if any of the source vectors in the initialized NMF model are redundant, 
#' and should be omitted from the actual NMF analysis
#' 
#' @param seed nmfModel object containing initialization of the factor matrices
#' @return boolean vector, indicating which source vector(s) are redundant
#' @importFrom NMF coef
#' @importFrom stats median
#' @author Nicolas Sauwen
checkForRedundantSources <- function(seed) {
  H0                   <- coef(seed)
  redundantSourceVect  <- rep(FALSE, nrow(H0))
  H0_sum               <- apply(H0,1,sum)
  redundantSourceVect[which(H0_sum < 0.005*median(H0_sum))] <- TRUE
  return(redundantSourceVect)
}

#' Remove redundant sources from the initial NMF model
#' 
#' @param seed nmfModel object containing initialization of the factor matrices
#' @param redundantSources boolean vector, obtained from \code{\link{checkForRedundantSources}}
#' @return nmfModel object with redundant sources removed from initial factor matrices
#' @importFrom NMF basis coef nmfModel
#' @author Nicolas Sauwen
removeRedundantSources <- function(seed, redundantSources) {
  W0 <- basis(seed)
  H0 <- coef(seed)
  redundantInds <- which(redundantSources == TRUE)
  W0 <- W0[,-redundantInds]
  H0 <- H0[-redundantInds,]
  seedNew <- nmfModel(W = W0 , H = H0)
  return(seedNew)
}

#' Re-introduce redundant source vectors and corresponding zero abundances 
#' into final NMF result
#' 
#' @param NMFResult Fitted NMF model
#' @param seed_orig Initial NMF model
#' @param redundantSources boolean vector, obtained from \code{\link{checkForRedundantSources}}
#' @return Final NMF model with redundant sources re-introduced
#' @importFrom NMF basis coef
#' @author Nicolas Sauwen
includeRedundantSources <- function(NMFResult, seed_orig, redundantSources) {
  W <- basis(NMFResult)
  H <- coef(NMFResult)
  W0 <- basis(seed_orig)
  H0 <- coef(seed_orig)
  nRows <- nrow(W0)
  nCols <- ncol(H0)
  rank <- ncol(W0)
  Wnew <- matrix(0, nrow = nRows, ncol = rank)
  Hnew <- matrix(0, nrow = rank, ncol = nCols)
  count <- 1
  for(i in 1:rank) {
    if(redundantSources[i] == TRUE){
      Wnew[,i] <- W0[,i]
      Hnew[i,] <- 0
    }
    else {
      Wnew[,i] <- W[,count]
      Hnew[i,] <- H[count,]
      count <- count + 1
    }
  }
  NMF::basis(NMFResult) <- Wnew
  NMF::coef(NMFResult) <- Hnew
  return(NMFResult)
}

#' Upsample NMF result to original temporal resolution
#' 
#' @param NMFResult Fitted NMF model
#' @param timePoints Original time points
#' @param subsamplingFactor Subsampling factor
#' @param shift Integer that correctly shifts subsampling index when 
#' applying NMF to multiple experiments
#' @return Upsampled NMF model
#' @importFrom stats approx
#' @author Nicolas Sauwen
upsampleNMFResult <- function(NMFResult, timePoints, subsamplingFactor, shift = 0) {
  H_subsamp <- coef(NMFResult)
  timePoints_subsamp <- timePoints[seq(1+shift,length(timePoints),subsamplingFactor)]
#  H <- apply(H_subsamp, 1, approx, x = timePoints_subsamp, xout = timePoints)
  H <- matrix(0, nrow(H_subsamp), length(timePoints))
  for(i in 1:nrow(H)) {
    H_interp <- approx(x = timePoints_subsamp, y = H_subsamp[i,], xout = timePoints, rule = 2)
    H[i,] <- H_interp$y
  }
  NMF::coef(NMFResult) <- H
  return(NMFResult)
}

#' complete source matrix
#' @importFrom stats runif
#' @keywords internal
completeSourceMatrix <- function(W0, rank) {
	
	extraComponents <- rank - ncol(W0)
	for(i in 1:extraComponents) {
		W0 <- cbind(W0, runif(nrow(W0)))
	}
	return(W0)
}

#' Based on previously obtained NMF result \code{NMFResult}, estimate coefficients for a new 
#' spectralData object \code{object} using non-negative least squares fitting. The result is 
#' returned as as an NMF model.
#' 
#' @param object \code{\link{SpectraInTime-class}}
#' @param NMFResult Fitted NMF model
#' @return Fitted non-negative least squares result in the form of an NMF model
#' @author nsauwen
#' @importFrom NMF basis
#' @importFrom nnls nnls
#' @export
predictNNLS <- function(object, NMFResult) {
	W <- NMF::basis(NMFResult)
	rank <- ncol(W)
	spectra <- getNMFInputMatrix( object )
	H <- matrix(0, rank, ncol(spectra))
	for(iSignal in 1:ncol(spectra)) {
		nlsFit <- nnls::nnls(W,spectra[,iSignal]) 
		H[,iSignal] <- stats::coefficients(nlsFit)
	}
	predictResult <- NMF::nmfModel(W = W, H = H)
}


#' Compute relative residual per observation of an NMF fit to a spectral data set
#' @param object \code{\link{SpectraInTime-class}}
#' @param NMFResult Fitted NMF model
#' @return Dataframe, containing time (observation) vector and residual vector
#' @author nsauwen
#' @importFrom hNMF residualNMF
#' @export
computeNMFResidu <- function(object, NMFResult) {
	
	timePts <- getTimePoints(object)
	spectra <- t(getSpectra(object))
	residu <- hNMF::residualNMF(spectra, NMFResult)
	
	return(data.frame(time=timePts, residu=residu))
	
}

Try the spectralAnalysis package in your browser

Any scripts or data that you put into this service are public.

spectralAnalysis documentation built on May 1, 2019, 10:25 p.m.