File: algorithms-siNMF.R

package info (click to toggle)
r-cran-nmf 0.23.0-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 3,344 kB
  • sloc: cpp: 680; ansic: 7; makefile: 2
file content (122 lines) | stat: -rw-r--r-- 3,311 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
# Implementation of siNMF from Badea (2008)
# 
# Author: Renaud Gaujoux
# Creation: 09 Jul 2012
###############################################################################

#' @include registry-algorithms.R
NULL

siNMF_R <- function(i, v, data, beta0=1, scale=TRUE, ...){
	
    # retrieve each factor
    w <- basis(data); h <- coef(data);
	# fixed terms
	nb <- nbterms(data); nc <- ncterms(data)
    
	if( i == 1 ){
	    if( !nc )
	  		stop("Method 'siNMF' requires a formula based model")
		
		if( !is.na(beta0) ){
			# compute beta
		    gr <- as.numeric(cterms(data)[[1L]])
		    beta <- beta0 * (norm(v[,gr==1], 'F') / norm(v[,gr==2], 'F'))^2
	    	# make sweeping vector
			vbeta <- rep(1, ncol(v))
			vbeta[gr==2] <- sqrt(beta)
		    staticVar('vbeta', vbeta, init=TRUE)
			# sweep data
			staticVar('v', sweep(v, 2L, vbeta, '*', check.margin=FALSE), init=TRUE)
		}
		# store non-fixed coef indexes
		staticVar('icoef', icoef(data), init=TRUE)
    }
    
    #precision threshold for numerical stability
    eps <- 10^-9
    
	sh <- h 
	if( !is.na(beta0) ){
		# retrieved swept matrix
		sv <- staticVar('v')
		vbeta <- staticVar('vbeta')
		# sweep h with beta
		sh <- sweep(h, 2L, vbeta, '*', check.margin=FALSE)
	}
	
    # compute standard euclidean updates
	w <- nmf_update.euclidean.w(sv, w, sh, eps=eps, nbterms=nb, ncterms=nc, copy=TRUE)
    h <- nmf_update.euclidean.h(v, w, h, eps=eps, nbterms=nb, ncterms=nc, copy=TRUE)
	
	# normalize columns of w
	if( scale ){
		icoef <- staticVar('icoef')
		wb <- w[, icoef]
		d <- sqrt(colSums(wb^2))
	    w[, icoef] <- sweep(wb, 2L, d, '/') 
		h[icoef, ] <- sweep(h[icoef, ], 1L, d, '*')
	}
	
	.coef(data) <- h
	.basis(data) <- w
    data
}

setNMFMethod('.siNMF', 'lee', Update=siNMF_R)

siNMF <- function(i, v, data, beta0=1, scale=TRUE, eps=10^-9, ...){
	
    # retrieve each factor
    w <- basis(data); h <- coef(data);
	# fixed terms
	nb <- nbterms(data); nc <- ncterms(data)
    
	if( i == 1 ){
	    if( !nc )
	  		stop("Method 'siNMF' requires a formula based model")
		
		vbeta <- NULL
		if( !is.na(beta0) ){
			# compute beta
		    gr <- cterms(data)[[1L]]
			gr <- droplevels(gr)
			# make sweeping vector
			vbeta <- rep(1, ncol(v))
			idx <- split(1:ncol(v), gr)
			# compute base value from first level
			beta <- beta0 * norm(v[,idx[[1]]], 'F')^2
			vbeta <- lapply(idx[-1], function(j){
				rep(beta / norm(v[,j], 'F')^2, length(j))
			})
			vbeta <- c(rep(1, length(idx[[1]])), unlist(vbeta, use.names=FALSE))
			vbeta <- vbeta[order(unlist(idx))]
		}
		# store weights
		staticVar('beta', vbeta, init=TRUE)
		# store non-fixed coef indexes
		staticVar('icoef', icoef(data), init=TRUE)
    }
    
    # retrieve weights
	beta <- staticVar('beta')
	
    # compute standard euclidean updates
	w <- nmf_update.euclidean.w(v, w, h, eps=eps, weight=beta, nbterms=nb, ncterms=nc, copy=FALSE)
    h <- nmf_update.euclidean.h(v, w, h, eps=eps, nbterms=nb, ncterms=nc, copy=FALSE)
	
	# normalize columns of w
	if( scale ){
		icoef <- staticVar('icoef')
		wb <- w[, icoef]
		d <- sqrt(colSums(wb^2))
	    w[, icoef] <- sweep(wb, 2L, d, '/', check.margin=FALSE) 
		h[icoef, ] <- sweep(h[icoef, ], 1L, d, '*', check.margin=FALSE)
	}
	
	.coef(data) <- h
	.basis(data) <- w
    data
}

nmfAlgorithm.siNMF <- setNMFMethod('siNMF', 'lee', Update=siNMF)