1#' @include NMFStrategy-class.R 2#' @include NMFfit-class.R 3NULL 4 5# Define union class for generalised function slots, e.g., slot 'NMFStrategyIterative::Stop' 6setClassUnion('.GfunctionSlotNULL', c('character', 'integer', 'numeric', 'function', 'NULL')) 7 8#' Interface for Algorithms: Implementation for Iterative NMF Algorithms 9#' 10#' @description 11#' This class provides a specific implementation for the generic function \code{run} 12#' -- concretising the virtual interface class \code{\linkS4class{NMFStrategy}}, 13#' for NMF algorithms that conform to the following iterative schema (starred numbers 14#' indicate mandatory steps): 15#' 16#' \itemize{ 17#' \item 1. Initialisation 18#' \item 2*. Update the model at each iteration 19#' \item 3. Stop if some criterion is satisfied 20#' \item 4. Wrap up 21#' } 22#' 23#' This schema could possibly apply to all NMF algorithms, since these are essentially optimisation algorithms, 24#' almost all of which use iterative methods to approximate a solution of the optimisation problem. 25#' The main advantage is that it allows to implement updates and stopping criterion separately, and combine them 26#' in different ways. 27#' In particular, many NMF algorithms are based on multiplicative updates, following the approach from 28#' \cite{Lee2001}, which are specially suitable to be cast into this simple schema. 29#' 30#' @slot onInit optional function that performs some initialisation or pre-processing on 31#' the model, before starting the iteration loop. 32#' @slot Update mandatory function that implement the update step, which computes new values for the model, based on its 33#' previous value. 34#' It is called at each iteration, until the stopping criterion is met or the maximum number of iteration is 35#' achieved. 36#' @slot Stop optional function that implements the stopping criterion. 37#' It is called \strong{before} each Update step. 38#' If not provided, the iterations are stopped after a fixed number of updates. 39#' @slot onReturn optional function that wraps up the result into an NMF object. 40#' It is called just before returning the 41#' 42setClass('NMFStrategyIterative' 43 , representation( 44 onInit = '.functionSlotNULL', 45 Update = '.functionSlot', # update method 46 Stop = '.GfunctionSlotNULL', # method called just after the update 47 onReturn = '.functionSlotNULL' # method called just before returning the resulting NMF object 48 ) 49 , prototype=prototype( 50 onInit = NULL 51 , Update = '' 52 , Stop = NULL 53 , onReturn = NULL 54 ) 55 , contains = 'NMFStrategy' 56 , validity = function(object){ 57 58 if( is.character(object@Update) && object@Update == '' ) 59 return("Slot 'Update' is required") 60 61 # check the arguments of methods 'Update' and 'Stop' 62 # (except for the 3 mandatory ones) 63 n.update <- names(formals(object@Update)) 64 65 # at least 3 arguments for 'Update' 66 if( length(n.update) < 3 ){ 67 return(str_c("Invalid 'Update' method - must have at least 3 arguments: ", 68 "current iteration number [i], ", 69 "target matrix [y], ", 70 "current NMF model iterate [x]")) 71 } 72 73 n.update <- n.update[-seq(3)] 74 # argument '...' must be present in method 'Update' 75 if( !is.element('...', n.update) ) 76 return("Invalid 'Update' method: must have argument '...' (even if not used)") 77 78 # at least 3 arguments for 'Stop' 79 if( !is.null(object@Stop) ){ 80 81 # retrieve the stopping criterion and check its intrinsic validity 82 .stop <- tryCatch( NMFStop(object@Stop, check=TRUE), 83 error = function(e) return(message(e))) 84 85 # Update and Stop methods cannot have overlapping arguments 86 n.stop <- names(formals(.stop)) 87 overlap <- intersect(n.update, n.stop) 88 overlap <- overlap[which(overlap!='...')] 89 if( length(overlap) > 0 ){ 90 return(str_c("Invalid 'Update' and 'Stop' methods: conflicting arguments ", 91 str_out(overlap, Inf))) 92 } 93 } 94 95 TRUE 96 } 97) 98 99 100#' Show method for objects of class \code{NMFStrategyIterative} 101#' @export 102setMethod('show', 'NMFStrategyIterative', 103 function(object){ 104 105 #cat('<object of class: NMFStrategyIterative>') 106 callNextMethod() 107 cat(" <Iterative schema>\n") 108 # go through the slots 109 s.list <- names(getSlots('NMFStrategyIterative')) 110 s.list <- setdiff(s.list, names(getSlots('NMFStrategy'))) 111 #s.list <- s.list[s.list=='ANY'] 112# s.list <- c('Update', 'Stop', 'WrapNMF') 113 out <- 114 sapply(s.list, function(sname){ 115 svalue <- slot(object,sname) 116 svalue <- 117 if( is.function(svalue) ) { 118 str_args(svalue, exdent=12) 119 } else if( is.null(svalue) ){ 120 'none' 121 } else { 122 paste("'", svalue,"'", sep='') 123 } 124 str_c(sname, ": ", svalue) 125 }) 126 cat(str_c(' ', out, collapse='\n'), "\n", sep='') 127 return(invisible()) 128 } 129) 130###% This class is an auxiliary class that defines the strategy's methods by directly callable functions. 131setClass('NMFStrategyIterativeX' 132 , contains = 'NMFStrategyIterative' 133 , representation = representation( 134 workspace = 'environment' # workspace to use persistent variables accross methods 135 ) 136) 137 138 139###% Creates a NMFStrategyIterativeX object from a NMFStrategyIterative object. 140xifyStrategy <- function(strategy, workspace=new.env(emptyenv())){ 141 142 # do nothing if already executable 143 if( is(strategy, 'NMFStrategyIterativeX') ) return(strategy) 144 145 # first check the strategy's validity 146 if( is.character(err <- validObject(strategy, test=TRUE)) ){ 147 stop("Invalid strategy definition:\n\t- ", err) 148 } 149 150 # intanciate the NMFStrategyIterativeX, creating the strategy's workspace 151 strategyX <- new('NMFStrategyIterativeX', strategy, workspace=workspace) 152 153 # define auxiliary function to preload the 'function' slots in class NMFStrategyIterativeX 154 preload.slot <- function(strategy, sname, default){ 155 156 # get the content of the slot 157 svalue <- slot(strategy,sname) 158 159 # if the slot is valid (i.e. it's a non-empty character string), then process the name into a valid function 160 fun <- 161 if( is.null(svalue) && !missing(default) ) default 162 else if( sname == 'Stop' ) NMFStop(svalue) 163 else if( is.character(svalue) && nchar(svalue) > 0 ){ 164 # set the slot with the executable version of the function 165 getFunction(svalue) 166 }else if( is.function(svalue) ) svalue 167 else 168 stop("NMFStrategyIterativeX - could not pre-load slot '", sname, "'") 169 170 # return the loaded function 171 fun 172 } 173 174 # preload the function slots 175 slot(strategyX, 'Update') <- preload.slot(strategyX, 'Update') 176 slot(strategyX, 'Stop') <- preload.slot(strategyX, 'Stop', function(strategy, i, target, data, ...){FALSE}) 177 slot(strategyX, 'onReturn') <- preload.slot(strategyX, 'onReturn', identity) 178 179 # load the objective function 180 objective(strategyX) <- nmfDistance(objective(strategy)) 181 182 # valid the preloaded object 183 validObject(strategyX) 184 185 # return the executable strategy 186 strategyX 187} 188 189# 190#setGeneric('Update', function(object, v, ...) standardGeneric('Update') ) 191#setMethod('Update', signature(object='NMFStrategyIterative', v='matrix'), function(object, v, ...){ object@data <- object@Update(v, object@data, ...) }) 192# 193#setGeneric('Stop', function(object, i) standardGeneric('Stop') ) 194#setMethod('Stop', signature(object='NMFStrategyIterative', i='integer'), function(object, i){ object@Stop(i, object@data) }) 195# 196#setGeneric('WrapNMF', function(object) standardGeneric('WrapNMF') ) 197#setMethod('WrapNMF', signature(object='NMFStrategyIterative'), function(object){ object@WrapNMF(object@data) }) 198 199###% Hook to initialize built-in iterative methods when the package is loaded 200 201 202###% Hook to initialize old R version built-in iterative methods 203 204#' Get/Set a Static Variable in NMF Algorithms 205#' 206#' @description 207#' This function is used in iterative NMF algorithms to manage variables 208#' stored in a local workspace, that are accessible to all functions that 209#' define the iterative schema described in \code{\linkS4class{NMFStrategyIterative}}. 210#' 211#' It is specially useful for computing stopping criteria, which often require model data from 212#' different iterations. 213#' 214#' @param name Name of the static variable (as a single character string) 215#' @param value New value of the static variable 216#' @param init a logical used when a \code{value} is provided, that specifies 217#' if the variable should be set to the new value only if it does not exist yet 218#' (\code{init=TRUE}). 219#' @return The value of the static variable 220#' @export 221staticVar <- local({ 222 223 .Workspace <- NULL 224 function(name, value, init=FALSE){ 225 226 # return last workspace 227 if( missing(name) ) return(.Workspace) 228 else if( is.null(name) ){ # reset workspace 229 .Workspace <<- NULL 230 return() 231 } else if( is.environment(name) ){ # setup up static environment 232 nmf.debug('Strategy Workspace', "initialize static workspace: ", capture.output(.Workspace), "=", capture.output(name)) 233 .Workspace <<- name 234 }else if( isString(name) && is.environment(.Workspace) ){ 235 if( missing(value) ){ 236 get(name, envir=.Workspace, inherits=FALSE) 237 }else{ 238 if( !init || !exists(name, envir=.Workspace, inherits=FALSE) ) 239 { 240 if( init ) nmf.debug('Strategy Workspace', "initialize variable '", name, "'") 241 assign(name, value, envir=.Workspace) 242 } 243 # return current value 244 get(name, envir=.Workspace, inherits=FALSE) 245 } 246 }else{ 247 stop("Invalid NMF workspace query: .Workspace=", class(.Workspace), '| name=', name 248 , if( !missing(value) ) paste0(' | value=', class(value))) 249 } 250 251 } 252}) 253 254#' Runs an NMF iterative algorithm on a target matrix \code{y}. 255#' 256#' @param .stop specification of a stopping criterion, that is used instead of the 257#' one associated to the NMF algorithm. 258#' It may be specified as: 259#' \itemize{ 260#' \item the access key of a registered stopping criterion; 261#' \item a single integer that specifies the exact number of iterations to perform, which will 262#' be honoured unless a lower value is explicitly passed in argument \code{maxIter}. 263#' \item a single numeric value that specifies the stationnarity threshold for the 264#' objective function, used in with \code{\link{nmf.stop.stationary}}; 265#' \item a function with signature \code{(object="NMFStrategy", i="integer", y="matrix", x="NMF", ...)}, 266#' where \code{object} is the \code{NMFStrategy} object that describes the algorithm being run, 267#' \code{i} is the current iteration, \code{y} is the target matrix and \code{x} is the current value of 268#' the NMF model. 269#' } 270#' @param maxIter maximum number of iterations to perform. 271#' 272#' @rdname NMFStrategy 273setMethod('run', signature(object='NMFStrategyIterative', y='matrix', x='NMFfit'), 274 function(object, y, x, .stop=NULL, maxIter = nmf.getOption('maxIter') %||% 2000L, ...){ 275 276 method <- object 277 # override the stop method on runtime 278 if( !is.null(.stop) ){ 279 method@Stop <- NMFStop(.stop) 280 # honour maxIter unless .stop is an integer and maxIter is not passed 281 # either directly or from initial call 282 # NB: maxIter may be not missing in the call to run() due to the application 283 # of default arguments from the Strategy within nmf(), in which case one does not 284 # want to honour it, since it is effectively missing in the original call. 285 if( is.integer(.stop) && (missing(maxIter) || !('maxIter' %in% names(x@call))) ) 286 maxIter <- .stop[1] 287 } 288 289 # debug object in debug mode 290 if( nmf.getOption('debug') ) show(method) 291 292 #Vc# Define local workspace for static variables 293 # this function can be called in the methods to get/set/initialize 294 # variables that are persistent within the strategy's workspace 295 .Workspace <- new.env() 296 staticVar(.Workspace) 297 on.exit( staticVar(NULL) ) 298 299 # runtime resolution of the strategy's functions by their names if necessary 300 strategyX = xifyStrategy(method, .Workspace) 301 run(strategyX, y, x, maxIter=maxIter, ...) 302}) 303 304#' @rdname NMFStrategy 305setMethod('run', signature(object='NMFStrategyIterativeX', y='matrix', x='NMFfit'), 306 function(object, y, x, maxIter, ...){ 307 308 strategy <- object 309 v <- y 310 seed <- x 311 #V!# NMFStrategyIterativeX::run 312 313 #Vc# Define workspace accessor function 314 # this function can be called in the methods to get/set/initialize 315 # variables that are persistent within the strategy's workspace 316# .Workspace <- strategy@workspace 317# assign('staticVar', function(name, value, init=FALSE){ 318# if( missing(value) ){ 319# get(name, envir=.Workspace, inherits=FALSE) 320# }else{ 321# if( !init || !exists(name, envir=.Workspace, inherits=FALSE) ) 322# { 323# if( init ) nmf.debug('Strategy Workspace', "initialize variable '", name, "'") 324# assign(name, value, envir=.Workspace) 325# } 326# } 327# } 328# , envir=.Workspace) 329 330 #Vc# initialize the strategy 331 # check validity of arguments if possible 332 method.args <- nmfFormals(strategy, runtime=TRUE) 333 internal.args <- method.args$internals 334 expected.args <- method.args$defaults 335 passed.args <- names(list(...)) 336 forbidden.args <- is.element(passed.args, c(internal.args)) 337 if( any(forbidden.args) ){ 338 stop("NMF::run - Update/Stop method : formal argument(s) " 339 , paste( paste("'", passed.args[forbidden.args],"'", sep=''), collapse=', ') 340 , " already set internally.", call.=FALSE) 341 } 342 # !is.element('...', expected.args) && 343 if( any(t <- !pmatch(passed.args, names(expected.args), nomatch=FALSE)) ){ 344 stop("NMF::run - onInit/Update/Stop method for algorithm '", name(strategy),"': unused argument(s) " 345 , paste( paste("'", passed.args[t],"'", sep=''), collapse=', '), call.=FALSE) 346 } 347 # check for required arguments 348 required.args <- sapply(expected.args, function(x){ x <- as.character(x); length(x) == 1 && nchar(x) == 0 } ) 349 required.args <- names(expected.args[required.args]) 350 required.args <- required.args[required.args!='...'] 351 352 if( any(t <- !pmatch(required.args, passed.args, nomatch=FALSE)) ) 353 stop("NMF::run - Update/Stop method for algorithm '", name(strategy),"': missing required argument(s) " 354 , paste( paste("'", required.args[t],"'", sep=''), collapse=', '), call.=FALSE) 355 356 #Vc# Start iterations 357 nmfData <- seed 358 # cache verbose level 359 verbose <- verbose(nmfData) 360 361 # clone the object to allow the updates to work in place 362 if( verbose > 1L ) 363 message("# Cloning NMF model seed ... ", appendLF=FALSE) 364 nmfFit <- clone(fit(nmfData)) 365 if( verbose > 1L ) 366 message("[", C.ptr(fit(nmfData)), " -> ", C.ptr(nmfFit), "]") 367 368 ## onInit 369 if( is.function(strategy@onInit) ){ 370 if( verbose > 1L ) message("# Step 1 - onInit ... ", appendLF=TRUE) 371 nmfFit <- strategy@onInit(strategy, v, nmfFit, ...) 372 if( verbose > 1L ) message("OK") 373 } 374 ## 375 376 # pre-load slots 377 updateFun <- strategy@Update 378 stopFun <- strategy@Stop 379 380 showNIter.step <- 50L 381 showNIter <- verbose && maxIter >= showNIter.step 382 if( showNIter ){ 383 ndIter <- nchar(as.character(maxIter)) 384 itMsg <- paste0('Iterations: %', ndIter, 'i', "/", maxIter) 385 cat(itMsgBck <- sprintf(itMsg, 0)) 386 itMsgBck <- nchar(itMsgBck) 387 } 388 i <- 0L 389 while( TRUE ){ 390 391 #Vc# Stopping criteria 392 # check convergence (generally do not stop for i=0L, but only initialise static variables 393 stop.signal <- stopFun(strategy, i, v, nmfFit, ...) 394 395 # if the strategy ask for stopping, then stop the iteration 396 if( stop.signal || i >= maxIter ) break; 397 398 # increment i 399 i <- i+1L 400 401 if( showNIter && (i==1L || i %% showNIter.step == 0L) ){ 402 cat(paste0(rep("\r", itMsgBck), sprintf(itMsg, i))) 403 } 404 405 #Vc# update the matrices 406 nmfFit <- updateFun(i, v, nmfFit, ...) 407 408 # every now and then track the error if required 409 nmfData <- trackError(nmfData, deviance(strategy, nmfFit, v, ...), niter=i) 410 411 } 412 if( showNIter ){ 413 ended <- if( stop.signal ) 'converged' else 'stopped' 414 cat("\nDONE (", ended, " at ",i,'/', maxIter," iterations)\n", sep='') 415 } 416 417 # force to compute last error if not already done 418 nmfData <- trackError(nmfData, deviance(strategy, nmfFit, v, ...), niter=i, force=TRUE) 419 420 # store the fitted model 421 fit(nmfData) <- nmfFit 422 423 #Vc# wrap up 424 # let the strategy build the result 425 nmfData <- strategy@onReturn(nmfData) 426 if( !inherits(nmfData, 'NMFfit') ){ 427 stop('NMFStrategyIterative[', name(strategy), ']::onReturn did not return a "NMF" instance [returned: "', class(nmfData), '"]') 428 } 429 430 # set the number of iterations performed 431 niter(nmfData) <- i 432 433 #return the result 434 nmf.debug('NMFStrategyIterativeX::run', 'Done') 435 invisible(nmfData) 436}) 437 438 439#' @export 440nmfFormals.NMFStrategyIterative <- function(x, runtime=FALSE, ...){ 441 442 strategy <- xifyStrategy(x) 443 # from run method 444 m <- getMethod('run', signature(object='NMFStrategyIterative', y='matrix', x='NMFfit')) 445 run.args <- allFormals(m)[-(1:3)] 446 # onInit 447 init.args <- if( is.function(strategy@onInit) ) formals(strategy@onInit) 448 # Update 449 update.args <- formals(strategy@Update) 450 # Stop 451 stop.args <- formals(strategy@Stop) 452 # spplit internals and 453 internal.args <- names(c(init.args[1:3], update.args[1:3], stop.args[1:4])) 454 expected.args <- c(init.args[-(1:3)], update.args[-(1:3)], stop.args[-(1:4)]) 455 456 if( runtime ){ 457 # prepend registered default arguments 458 expected.args <- expand_list(strategy@defaults, expected.args) 459 list(internal=internal.args, defaults=expected.args) 460 }else{ 461 args <- c(run.args, expected.args) 462 # prepend registered default arguments 463 expand_list(strategy@defaults, args) 464 } 465} 466 467################################################################################################ 468# INITIALIZATION METHODS 469################################################################################################ 470 471################################################################################################ 472# UPDATE METHODS 473################################################################################################ 474 475#' NMF Multiplicative Updates for Kullback-Leibler Divergence 476#' 477#' Multiplicative updates from \cite{Lee2001} for standard Nonnegative Matrix Factorization 478#' models \eqn{V \approx W H}, where the distance between the target matrix and its NMF 479#' estimate is measured by the Kullback-Leibler divergence. 480#' 481#' \code{nmf_update.KL.w} and \code{nmf_update.KL.h} compute the updated basis and coefficient 482#' matrices respectively. 483#' They use a \emph{C++} implementation which is optimised for speed and memory usage. 484#' 485#' @details 486#' The coefficient matrix (\code{H}) is updated as follows: 487#' \deqn{ 488#' H_{kj} \leftarrow H_{kj} \frac{\left( sum_i \frac{W_{ik} V_{ij}}{(WH)_{ij}} \right)}{ sum_i W_{ik} }. 489#' }{ 490#' H_kj <- H_kj ( sum_i [ W_ik V_ij / (WH)_ij ] ) / ( sum_i W_ik ) 491#' } 492#' 493#' These updates are used in built-in NMF algorithms \code{\link[=KL-nmf]{KL}} and 494#' \code{\link[=brunet-nmf]{brunet}}. 495#' 496#' @param v target matrix 497#' @param w current basis matrix 498#' @param h current coefficient matrix 499#' @param nbterms number of fixed basis terms 500#' @param ncterms number of fixed coefficient terms 501#' @param copy logical that indicates if the update should be made on the original 502#' matrix directly (\code{FALSE}) or on a copy (\code{TRUE} - default). 503#' With \code{copy=FALSE} the memory footprint is very small, and some speed-up may be 504#' achieved in the case of big matrices. 505#' However, greater care should be taken due the side effect. 506#' We recommend that only experienced users use \code{copy=TRUE}. 507#' 508#' @return a matrix of the same dimension as the input matrix to update 509#' (i.e. \code{w} or \code{h}). 510#' If \code{copy=FALSE}, the returned matrix uses the same memory as the input object. 511#' 512#' @author 513#' Update definitions by \cite{Lee2001}. 514#' 515#' C++ optimised implementation by Renaud Gaujoux. 516#' 517#' @rdname nmf_update_KL 518#' @aliases nmf_update.KL 519#' @export 520nmf_update.KL.h <- std.divergence.update.h <- function(v, w, h, nbterms=0L, ncterms=0L, copy=TRUE) 521{ 522 .Call("divergence_update_H", v, w, h, nbterms, ncterms, copy, PACKAGE='NMF') 523} 524#' \code{nmf_update.KL.w_R} and \code{nmf_update.KL.h_R} implement the same updates 525#' in \emph{plain R}. 526#' 527#' @param wh already computed NMF estimate used to compute the denominator term. 528#' 529#' @rdname nmf_update_KL 530#' @export 531nmf_update.KL.h_R <- R_std.divergence.update.h <- function(v, w, h, wh=NULL) 532{ 533 # compute WH if necessary 534 if( is.null(wh) ) wh <- w %*% h 535 536 # divergence-reducing NMF iterations 537 # H_au = H_au ( sum_i [ W_ia V_iu / (WH)_iu ] ) / ( sum_k W_ka ) -> each row of H is divided by a the corresponding colSum of W 538 h * crossprod(w, v / wh) / colSums(w) 539} 540 541#' @details 542#' The basis matrix (\code{W}) is updated as follows: 543#' \deqn{ 544#' W_{ik} \leftarrow W_{ik} \frac{ sum_j [\frac{H_{kj} A_{ij}}{(WH)_{ij}} ] }{sum_j H_{kj} } 545#' }{ 546#' W_ik <- W_ik ( sum_u [H_kl A_il / (WH)_il ] ) / ( sum_l H_kl ) 547#' } 548#' @rdname nmf_update_KL 549#' @export 550nmf_update.KL.w <- std.divergence.update.w <- function(v, w, h, nbterms=0L, ncterms=0L, copy=TRUE) 551{ 552 .Call("divergence_update_W", v, w, h, nbterms, ncterms, copy, PACKAGE='NMF') 553} 554#' @rdname nmf_update_KL 555#' @export 556nmf_update.KL.w_R <- R_std.divergence.update.w <- function(v, w, h, wh=NULL) 557{ 558 # compute WH if necessary 559 if( is.null(wh) ) wh <- w %*% h 560 561 # W_ia = W_ia ( sum_u [H_au A_iu / (WH)_iu ] ) / ( sum_v H_av ) -> each column of W is divided by a the corresponding rowSum of H 562 #x2 <- matrix(rep(rowSums(h), nrow(w)), ncol=ncol(w), byrow=TRUE); 563 #w * tcrossprod(v / wh, h) / x2; 564 sweep(w * tcrossprod(v / wh, h), 2L, rowSums(h), "/", check.margin = FALSE) # optimize version? 565 566} 567 568 569 570#' NMF Multiplicative Updates for Euclidean Distance 571#' 572#' Multiplicative updates from \cite{Lee2001} for standard Nonnegative Matrix Factorization 573#' models \eqn{V \approx W H}, where the distance between the target matrix and its NMF 574#' estimate is measured by the -- euclidean -- Frobenius norm. 575#' 576#' \code{nmf_update.euclidean.w} and \code{nmf_update.euclidean.h} compute the updated basis and coefficient 577#' matrices respectively. 578#' They use a \emph{C++} implementation which is optimised for speed and memory usage. 579#' 580#' @details 581#' The coefficient matrix (\code{H}) is updated as follows: 582#' \deqn{ 583#' H_{kj} \leftarrow \frac{\max(H_{kj} W^T V)_{kj}, \varepsilon) }{(W^T W H)_{kj} + \varepsilon} 584#' }{ 585#' H_kj <- max(H_kj (W^T V)_kj, eps) / ( (W^T W H)_kj + eps ) 586#' } 587#' 588#' These updates are used by the built-in NMF algorithms \code{\link[=Frobenius-nmf]{Frobenius}} and 589#' \code{\link[=lee-nmf]{lee}}. 590#' 591#' @inheritParams nmf_update.KL.h 592#' @param eps small numeric value used to ensure numeric stability, by shifting up 593#' entries from zero to this fixed value. 594#' 595#' @return a matrix of the same dimension as the input matrix to update 596#' (i.e. \code{w} or \code{h}). 597#' If \code{copy=FALSE}, the returned matrix uses the same memory as the input object. 598#' 599#' @author 600#' Update definitions by \cite{Lee2001}. 601#' 602#' C++ optimised implementation by Renaud Gaujoux. 603#' 604#' @rdname nmf_update_euclidean 605#' @aliases nmf_update.euclidean 606#' @export 607nmf_update.euclidean.h <- std.euclidean.update.h <- 608function(v, w, h, eps=10^-9, nbterms=0L, ncterms=0L, copy=TRUE){ 609 .Call("euclidean_update_H", v, w, h, eps, nbterms, ncterms, copy, PACKAGE='NMF') 610} 611#' \code{nmf_update.euclidean.w_R} and \code{nmf_update.euclidean.h_R} implement the same updates 612#' in \emph{plain R}. 613#' 614#' @param wh already computed NMF estimate used to compute the denominator term. 615#' 616#' @rdname nmf_update_euclidean 617#' @export 618nmf_update.euclidean.h_R <- R_std.euclidean.update.h <- function(v, w, h, wh=NULL, eps=10^-9){ 619 # compute WH if necessary 620 den <- if( is.null(wh) ) crossprod(w) %*% h 621 else{ t(w) %*% wh} 622 623 # H_au = H_au (W^T V)_au / (W^T W H)_au 624 pmax(h * crossprod(w,v),eps) / (den + eps); 625} 626 627#' @details 628#' The basis matrix (\code{W}) is updated as follows: 629#' \deqn{ 630#' W_ik \leftarrow \frac{\max(W_ik (V H^T)_ik, \varepsilon) }{ (W H H^T)_ik + \varepsilon} 631#' }{ 632#' W_ik <- max(W_ik (V H^T)_ik, eps) / ( (W H H^T)_ik + eps ) 633#' } 634#' 635#' @param weight numeric vector of sample weights, e.g., used to normalise samples 636#' coming from multiple datasets. 637#' It must be of the same length as the number of samples/columns in \code{v} 638#' -- and \code{h}. 639#' 640#' @rdname nmf_update_euclidean 641#' @export 642nmf_update.euclidean.w <- std.euclidean.update.w <- 643function(v, w, h, eps=10^-9, nbterms=0L, ncterms=0L, weight=NULL, copy=TRUE){ 644 .Call("euclidean_update_W", v, w, h, eps, weight, nbterms, ncterms, copy, PACKAGE='NMF') 645} 646#' @rdname nmf_update_euclidean 647#' @export 648nmf_update.euclidean.w_R <- R_std.euclidean.update.w <- function(v, w, h, wh=NULL, eps=10^-9){ 649 # compute WH if necessary 650 den <- if( is.null(wh) ) w %*% tcrossprod(h) 651 else{ wh %*% t(h)} 652 653 # W_ia = W_ia (V H^T)_ia / (W H H^T)_ia and columns are rescaled after each iteration 654 pmax(w * tcrossprod(v, h), eps) / (den + eps); 655} 656 657 658################################################################################################ 659# AFTER-UPDATE METHODS 660################################################################################################ 661 662#' Stopping Criteria for NMF Iterative Strategies 663#' 664#' The function documented here implement stopping/convergence criteria 665#' commonly used in NMF algorithms. 666#' 667#' \code{NMFStop} acts as a factory method that creates stopping criterion functions 668#' from different types of values, which are subsequently used by 669#' \code{\linkS4class{NMFStrategyIterative}} objects to determine when to stop their 670#' iterative process. 671#' 672#' @details 673#' \code{NMFStop} can take the following values: 674#' \describe{ 675#' \item{function}{ is returned unchanged, except when it has no arguments, 676#' in which case it assumed to be a generator, which is immediately called and should return 677#' a function that implements the actual stopping criterion;} 678#' \item{integer}{ the value is used to create a stopping criterion that stops at 679#' that exact number of iterations via \code{nmf.stop.iteration};} 680#' \item{numeric}{ the value is used to create a stopping criterion that stops when 681#' at that stationary threshold via \code{nmf.stop.threshold};} 682#' \item{character}{ must be a single string which must be an access key 683#' for registered criteria (currently available: \dQuote{connectivity} and \dQuote{stationary}), 684#' or the name of a function in the global environment or the namespace of the loading package.} 685#' } 686#' 687#' @param s specification of the stopping criterion. 688#' See section \emph{Details} for the supported formats and how they are processed. 689#' @param check logical that indicates if the validity of the stopping criterion 690#' function should be checked before returning it. 691#' 692#' @return a function that can be passed to argument \code{.stop} of function 693#' \code{\link{nmf}}, which is typically used when the algorith is implemented as 694#' an iterative strategy. 695#' 696#' @aliases stop-NMF 697#' @rdname stop-NMF 698#' @export 699NMFStop <- function(s, check=TRUE){ 700 701 key <- s 702 703 fun <- 704 if( is.integer(key) ) nmf.stop.iteration(key) 705 else if( is.numeric(key) ) nmf.stop.threshold(key) 706 else if( is.function(key) ) key 707 else if( is.character(key) ){ 708 # update .stop for back compatibility: 709 if( key == 'nmf.stop.consensus') key <- 'connectivity' 710 711 # first lookup for a `nmf.stop.*` function 712 key2 <- paste('nmf.stop.', key, sep='') 713 e <- pkgmaker::packageEnv() 714 sfun <- getFunction(key2, mustFind=FALSE, where = e) 715 if( is.null(sfun) ) # lookup for the function as such 716 sfun <- getFunction(key, mustFind = FALSE, where = e) 717 if( is.null(sfun) ) 718 stop("Invalid key ['", key,"']: could not find functions '",key2, "' or '", key, "'") 719 sfun 720 }else if( identical(key, FALSE) ) # create a function that does not stop 721 function(strategy, i, target, data, ...){FALSE} 722 else 723 stop("Invalid key: should be a function, a character string or a single integer/numeric value. See ?NMFStop.") 724 725 # execute if generator (i.e. no arguments) 726 if( length(formals(fun)) == 0L ) fun <- fun() 727 728 # check validity if requested 729 if( check ){ 730 n.stop <- names(formals(fun)) 731 if( length(n.stop) < 4 ){ 732 stop("Invalid 'Stop' method - must have at least 4 arguments: ", 733 "NMF strategy object [strategy], ", 734 "current iteration number [i], ", 735 "target matrix [y], ", 736 "current NMF model iterate [x]") 737 } 738 739 n.stop <- n.stop[-seq(4)] 740 # argument '...' must be present in method 'Stop' 741 if( !is.element('...', n.stop) ) 742 stop("Invalid 'Stop' method: must have argument '...' (even if not used)") 743 } 744 745 # return function 746 fun 747} 748 749#' \code{nmf.stop.iteration} generates a function that implements the stopping 750#' criterion that limits the number of iterations to a maximum of \code{n}), 751#' i.e. that returns \code{TRUE} if \code{i>=n}, \code{FALSE} otherwise. 752#' 753#' @param n maximum number of iteration to perform. 754#' 755#' @return a function that can be used as a stopping criterion for NMF algorithms 756#' defined as \code{\linkS4class{NMFStrategyIterative}} objects. 757#' That is a function with arguments \code{(strategy, i, target, data, ...)} 758#' that returns \code{TRUE} if the stopping criterion is satisfied -- which in 759#' turn stops the iterative process, and \code{FALSE} otherwise. 760#' 761#' @export 762#' @family NMFStrategyIterative 763#' @rdname stop-NMF 764nmf.stop.iteration <- function(n){ 765 766 nmf.debug("Using stopping criterion - Fixed number of iterations: ", n) 767 if( !is.numeric(n) ) 768 stop("Invalid argument `n`: must be an integer value") 769 if( length(n) > 1 ) 770 warning("NMF::nmf - Argument `n` [", deparse(substitute(n)), "] has length > 1: only using the first element.") 771 772 .max <- n[1] 773 function(object, i, y, x, ...) i >= .max 774} 775 776#' \code{nmf.stop.threshold} generates a function that implements the stopping 777#' criterion that stops when a given stationarity threshold is achieved by 778#' successive iterations. 779#' The returned function is identical to \code{nmf.stop.stationary}, but with 780#' the default threshold set to \code{threshold}. 781#' 782#' @param threshold default stationarity threshold 783#' 784#' @export 785#' @rdname stop-NMF 786nmf.stop.threshold <- function(threshold){ 787 788 nmf.debug("Using stopping criterion - Stationarity threshold: ", threshold) 789 if( !is.numeric(threshold) ) 790 stop("Invalid argument `threshold`: must be a numeric value") 791 if( length(threshold) > 1 ) 792 warning("NMF::nmf - Argument `threshold` [", deparse(substitute(threshold)), "] has length > 1: only using the first element.") 793 794 eval(parse(text=paste("function(strategy, i, target, data, stationary.th=", threshold, ", ...){ 795 nmf.stop.stationary(strategy, i, target, data, stationary.th=stationary.th, ...) 796 }"))) 797} 798 799 800#' \code{nmf.stop.stationary} implements the stopping criterion of stationarity 801#' of the objective value, which stops when the gradient of the objective function 802#' is uniformly small over a certain number of iterations. 803#' 804#' More precisely, the objective function is computed over \eqn{n} successive iterations (specified 805#' in argument \code{check.niter}), every \code{check.interval} iterations. 806#' The criterion stops when the absolute difference between the maximum and the minimum 807#' objective values over these iterations is lower than a given threshold \eqn{\alpha} 808#' (specified in \code{stationary.th}): 809#' 810#' \deqn{ 811#' \left| \frac{\max_{i- N_s + 1 \leq k \leq i} D_k - \min_{i - N_s +1 \leq k \leq i} D_k}{n} \right| \leq \alpha, 812#' }{ 813#' | [max( D(i- N_s + 1), ..., D(i) ) - min( D(i- N_s + 1), ..., D(i) )] / n | <= alpha 814#' } 815#' 816#' @param object an NMF strategy object 817#' @param i the current iteration 818#' @param y the target matrix 819#' @param x the current NMF model 820#' @param stationary.th maximum absolute value of the gradient, for the objective 821#' function to be considered stationary. 822#' @param check.interval interval (in number of iterations) on which the stopping 823#' criterion is computed. 824#' @param check.niter number of successive iteration used to compute the stationnary 825#' criterion. 826#' @param ... extra arguments passed to the function \code{\link{objective}}, 827#' which computes the objective value between \code{x} and \code{y}. 828#' 829#' @export 830#' @rdname stop-NMF 831nmf.stop.stationary <- local({ 832 833 # static variable 834 .last.objective.value <- c(-Inf, Inf) 835 .niter <- 0L 836 837 .store_value <- function(value){ 838 .niter <<- .niter + 1L 839 .last.objective.value <<- c(max(.last.objective.value[1L], value) 840 , min(.last.objective.value[2L], value)) 841 } 842 843 .reset_value <- function(){ 844 .last.objective.value <<- c(-Inf, Inf) 845 .niter <<- 0L 846 } 847 848 function(object, i, y, x, stationary.th=.Machine$double.eps, check.interval=5*check.niter, check.niter=10L, ...){ 849 850 # check validity 851 if( check.interval < check.niter ){ 852 stop("Invalid argument values: `check.interval` must always be greater than `check.niter`") 853 } 854 # initialisation call: compute initial objective value 855 if( i == 0L || (i == 1L && is.null(.last.objective.value)) ){ 856 .reset_value() 857 858 # give the chance to update once and estimate from a partial model 859 if( is.partial.nmf(x) ) return( FALSE ) 860 861 # compute initial deviance 862 current.value <- deviance(object, x, y, ...) 863 # check for NaN, i.e. probably infinitely small value (cf. bug reported by Nadine POUKEN SIEWE) 864 if( is.nan(current.value) ) return(TRUE) 865 866 # store value in static variable for next calls 867 .store_value(current.value) 868 869 return(FALSE) 870 } 871 872 # test convergence only every 10 iterations 873 if( .niter==0L && i %% check.interval != 0 ) return( FALSE ); 874 875 # get last objective value from static variable 876 current.value <- deviance(object, x, y, ...) 877 # check for NaN, i.e. probably infinitely small value (cf. bug reported by Nadine POUKEN SIEWE) 878 if( is.nan(current.value) ) return(TRUE) 879 880 # update static variables 881 .store_value(current.value) 882 883 # once values have been computed for check.niter iterations: 884 # check if the difference in the extreme objective values is small enough 885 if( .niter == check.niter+1 ){ 886 crit <- abs(.last.objective.value[1L] - .last.objective.value[2L]) / check.niter 887 if( crit <= stationary.th ){ 888 if( nmf.getOption('verbose') ){ 889 message(crit) 890 } 891 return( TRUE ) 892 } 893 .reset_value() 894 } 895 896 # do NOT stop 897 FALSE 898 } 899}) 900 901#' \code{nmf.stop.connectivity} implements the stopping criterion that is based 902#' on the stationarity of the connectivity matrix. 903#' 904#' @inheritParams nmf.stop.stationary 905#' @param stopconv number of iterations intervals over which the connectivity 906#' matrix must not change for stationarity to be achieved. 907#' 908#' @export 909#' @rdname stop-NMF 910nmf.stop.connectivity <- local({ 911 912 # static variables 913 .consold <- NULL 914 .inc <- NULL 915 916 function(object, i, y, x, stopconv=40, check.interval=10, ...){ 917 918 if( i == 0L ){ # initialisation call 919 # Initialize consensus variables 920 # => they are static variables within the strategy's workspace so that 921 # they are persistent and available throughout across the calls 922 p <- ncol(x) 923 .consold <<- matrix(0, p, p) 924 .inc <<- 0 925 return(FALSE) 926 } 927 928 # test convergence only every 10 iterations 929 if( i %% check.interval != 0 ) return( FALSE ); 930 931 # retrieve metaprofiles 932 h <- coef(x, all=FALSE) 933 934 # construct connectivity matrix 935 index <- apply(h, 2, function(x) which.max(x) ) 936 cons <- outer(index, index, function(x,y) ifelse(x==y, 1,0)); 937 938 changes <- cons != .consold 939 if( !any(changes) ) .inc <<- .inc + 1 # connectivity matrix has not changed: increment the count 940 else{ 941 .consold <<- cons; 942 .inc <<- 0; # else restart counting 943 } 944 945 # prints number of changing elements 946 #if( verbose(x) ) cat( sprintf('%d ', sum(changes)) ) 947 #cat( sprintf('%d ', sum(changes)) ) 948 949 # assume convergence is connectivity stops changing 950 if( .inc > stopconv ) return( TRUE ); 951 952 # do NOT stop 953 FALSE 954 } 955}) 956 957 958################################################################################################ 959# WRAP-UP METHODS 960################################################################################################ 961