Diff of /pkg/R/Matrix.R

Parent Directory | Revision Log | Patch

	revision 2110, Sat Jan 26 20:59:26 2008 UTC	revision 2115, Sat Feb 23 09:23:17 2008 UTC
#	Line 45	Line 45
45	function(x, ...) as.logical(as.vector(x)))	function(x, ...) as.logical(as.vector(x)))
46
47	setMethod("cov2cor", signature(V = "Matrix"),	setMethod("cov2cor", signature(V = "Matrix"),
48	function(V) as(cov2cor(as(V, "matrix")), "dpoMatrix"))	function(V) { ## was as(cov2cor(as(V, "matrix")), "dpoMatrix"))
49		r <- V
50		p <- (d <- dim(V))[1]
51		if(p != d[2]) stop("'V' is not a square matrix")
52		Is <- sqrt(1/diag(V)) # diag( 1/sigma_i )
53		if(any(!is.finite(Is)))
54		warning("diag(.) had 0 or NA entries; non-finite result is doubtful")
55		Is <- Diagonal(x = Is)
56		r <- Is %*% V %*% Is
57		r[cbind(1L:p,1L:p)] <- 1 # exact in diagonal
58		as(forceSymmetric(r), "dpoMatrix")
59		})
60
61	## "base" has an isSymmetric() S3-generic since R 2.3.0	## "base" has an isSymmetric() S3-generic since R 2.3.0
62	setMethod("isSymmetric", signature(object = "symmetricMatrix"),	setMethod("isSymmetric", signature(object = "symmetricMatrix"),
63	function(object,tol) TRUE)	function(object, ...) TRUE)
64	setMethod("isSymmetric", signature(object = "triangularMatrix"),	setMethod("isSymmetric", signature(object = "triangularMatrix"),
65	## TRUE iff diagonal:	## TRUE iff diagonal:
66	function(object,tol) isDiagonal(object))	function(object, ...) isDiagonal(object))
	setMethod("isTriangular", signature(object = "triangularMatrix"),
	function(object, ...) TRUE)
67
68	setMethod("isTriangular", signature(object = "matrix"), isTriMat)	setMethod("isTriangular", signature(object = "matrix"), isTriMat)
69
70	setMethod("isDiagonal", signature(object = "matrix"), .is.diagonal)	setMethod("isDiagonal", signature(object = "matrix"), .is.diagonal)
71
72		## The "catch all" methods -- far from optimal:
73		setMethod("symmpart", signature(x = "Matrix"),
74		function(x) as((x + t(x))/2, "symmetricMatrix"))
75		setMethod("skewpart", signature(x = "Matrix"),
76		function(x) (x - t(x))/2)
77
78		## FIXME: do this (similarly as for "ddense.." in C
79		setMethod("symmpart", signature(x = "matrix"), function(x) (x + t(x))/2)
80		setMethod("skewpart", signature(x = "matrix"), function(x) (x - t(x))/2)
81
82
83
84
85	setMethod("dim", signature(x = "Matrix"),	setMethod("dim", signature(x = "Matrix"),
#	Line 101	Line 120
120	##        "!" is in ./not.R	##        "!" is in ./not.R
121
122
123	Matrix <-	Matrix <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE,
124	function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL,	dimnames = NULL, sparse = NULL, forceCheck = FALSE)
	sparse = NULL, forceCheck = FALSE)
125	{	{
126	sparseDefault <- function(m) prod(dim(m)) > 2*sum(isN0(as(m, "matrix")))	sparseDefault <- function(m) prod(dim(m)) > 2*sum(isN0(as(m, "matrix")))
127
#	Line 112	Line 130
130	class(data) <- "matrix" # "matrix" first for S4 dispatch	class(data) <- "matrix" # "matrix" first for S4 dispatch
131	if(is.null(sparse1 <- sparse) && (i.M || is(data, "matrix")))	if(is.null(sparse1 <- sparse) && (i.M || is(data, "matrix")))
132	sparse <- sparseDefault(data)	sparse <- sparseDefault(data)
133		sM <- FALSE
134	doDN <- TRUE	doDN <- TRUE
135	if (i.M) {	if (i.M) {
136	if(!missing(nrow) || !missing(ncol)|| !missing(byrow))	if(!missing(nrow) || !missing(ncol)|| !missing(byrow))
#	Line 131	Line 149
149	## Matrix(0, ...) : always sparse unless "sparse = FALSE":	## Matrix(0, ...) : always sparse unless "sparse = FALSE":
150	if(is.null(sparse)) sparse1 <- sparse <- TRUE	if(is.null(sparse)) sparse1 <- sparse <- TRUE
151	i.M <- sM <- TRUE	i.M <- sM <- TRUE
152		isSym <- nrow == ncol
153	## will be sparse: do NOT construct full matrix!	## will be sparse: do NOT construct full matrix!
154	data <- new(if(is.numeric(data)) "dgTMatrix" else	data <- new(paste(if(is.numeric(data)) "d" else
155	if(is.logical(data)) "lgTMatrix" else	if(is.logical(data)) "l" else
156	stop("invalid 'data'"),	stop("invalid 'data'"),
157		if(isSym) "s" else "g", "CMatrix", sep=''),
158		p = rep.int(0L, ncol+1L),
159	Dim = as.integer(c(nrow,ncol)),	Dim = as.integer(c(nrow,ncol)),
160	Dimnames = if(is.null(dimnames)) list(NULL,NULL)	Dimnames = if(is.null(dimnames)) list(NULL,NULL)
161	else dimnames)	else dimnames)
#	Line 162	Line 183
183	isTri <- isTriangular(data)	isTri <- isTriangular(data)
184	isDiag <- isSym # cannot be diagonal if it isn't symmetric	isDiag <- isSym # cannot be diagonal if it isn't symmetric
185	if(isDiag)	if(isDiag)
186	isDiag <- isDiagonal(data)	isDiag <- !isTRUE(sparse1) && isDiagonal(data)
187
188	## Find proper matrix class 'cl'	## try to coerce ``via'' virtual classes
189	cl <-	if(isDiag) { ## diagonal is preferred to sparse !
190	if(isDiag && !isTRUE(sparse1))	data <- as(data, "diagonalMatrix")
191	"diagonalMatrix" # -> will automatically check for type	isSym <- FALSE
192	else {	} else if(sparse && !sM)
	## consider it's type
	ctype <-
	if(is(data,"Matrix")) class(data)
	else {
	if("complex" == (ctype <- typeof(data)))
	"z" else ctype
	}
	ctype <- substr(ctype, 1,1) # "d", "l", "i" or "z"
	if(ctype == "z")
	stop("complex matrices not yet implemented in Matrix package")
	if(ctype == "i") {
	warning("integer matrices not yet implemented in 'Matrix'; ",
	"using 'double' ones'")
	ctype <- "d"
	}
	paste(ctype,
	if(sparse) {
	if(isSym) "sCMatrix" else
	if(isTri) "tCMatrix" else "gCMatrix"
	} else { ## dense
	if(isSym) "syMatrix" else
	if(isTri) "trMatrix" else "geMatrix"
	}, sep="")
	}
	## Can we coerce and be done?
	if(!canCoerce(data,cl)) { ## try to coerce ``via'' virtual classes
	if(sparse && !sM)
193	data <- as(data, "sparseMatrix")	data <- as(data, "sparseMatrix")
194	else if(!sparse && !is(data, "denseMatrix"))	else if(!sparse) {
195		if(i.M) { ## data is 'Matrix'
196		if(!is(data, "denseMatrix"))
197	data <- as(data, "denseMatrix")	data <- as(data, "denseMatrix")
198	if(isTri && !is(data, "triangularMatrix"))	} else { ## data is "matrix" (and result "dense" -> go via "general"
199	data <- as(data, "triangularMatrix")	ctype <- typeof(data)
200	else if(isSym && !is(data, "symmetricMatrix"))	if (ctype == "complex")
201	data <- as(data, "symmetricMatrix")	stop("complex matrices not yet implemented in Matrix package")
202		if (ctype == "integer") ## integer Matrices not yet implemented
203		storage.mode(data) <- "double"
204		data <- new(paste(.M.kind(data), "geMatrix", sep=''),
205		Dim = dim(data),
206		Dimnames = .M.DN(data),
207		x = c(data))
208	}	}
209	## now coerce in any case .. maybe producing sensible error message:	}
210	as(data, cl)
211		if(isTri && !is(data, "triangularMatrix")) {
212		data <- if(attr(isTri,"kind") == "L") tril(data) else triu(data)
213		#was as(data, "triangularMatrix")
214		} else if(isSym && !is(data, "symmetricMatrix"))
215		data <- forceSymmetric(data) #was as(data, "symmetricMatrix")
216
217		data
218	}	}
219
220	## Methods for operations where one argument is numeric	## Methods for operations where one argument is numeric
#	Line 376	Line 384
384
385	## missing 'drop' --> 'drop = TRUE'	## missing 'drop' --> 'drop = TRUE'
386	##                     -----------	##                     -----------
387	## select rows	## select rows __ or __ vector indexing:
388	setMethod("[", signature(x = "Matrix", i = "index", j = "missing",	setMethod("[", signature(x = "Matrix", i = "index", j = "missing",
389	drop = "missing"),	drop = "missing"),
390	function(x,i,j, ..., drop) {	function(x,i,j, ..., drop) {
391	if(nargs() == 2) { ## e.g. M[0] , M[TRUE],  M[1:2]	if(nargs() == 2) { ## e.g. M[0] , M[TRUE],  M[1:2]
392	if(any(i) || prod(dim(x)) == 0)	if(any(as.logical(i)) || prod(dim(x)) == 0)
393		## FIXME: for *large sparse*, use sparseVector !
394	as.vector(x)[i]	as.vector(x)[i]
395	else ## save memory	else ## save memory (for large sparse M):
396	as.vector(x[1,1])[FALSE]	as.vector(x[1,1])[FALSE]
397	} else {	} else {
398	callGeneric(x, i=i, , drop=TRUE)	callGeneric(x, i=i, , drop=TRUE)
#	Line 412	Line 421
421	nA <- nargs()	nA <- nargs()
422	if(nA == 2) { ##  M [ M >= 7 ]	if(nA == 2) { ##  M [ M >= 7 ]
423	## FIXME: when both 'x' and 'i' are sparse, this can be very inefficient	## FIXME: when both 'x' and 'i' are sparse, this can be very inefficient
424	as(x, geClass(x))@x[as.vector(i)]	if(is(x, "sparseMatrix"))
425		message("<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient")
426		toC <- geClass(x)
427		if(canCoerce(x, toC)) as(x, toC)@x[as.vector(i)]
428		else as(as(as(x, "generalMatrix"), "denseMatrix"), toC)@x[as.vector(i)]
429	## -> error when lengths don't match	## -> error when lengths don't match
430	} else if(nA == 3) { ##  M [ M[,1, drop=FALSE] >= 7, ]	} else if(nA == 3) { ##  M [ M[,1, drop=FALSE] >= 7, ]
431	stop("not-yet-implemented 'Matrix' subsetting") ## FIXME	stop("not-yet-implemented 'Matrix' subsetting") ## FIXME
#	Line 428	Line 441
441	.M.sub.i.logical)	.M.sub.i.logical)
442
443
444		subset.ij <- function(x, ij) {
445		m <- nrow(ij)
446		if(m > 3) {
447		cld <- getClassDef(class(x))
448		sym.x <- extends(cld, "symmetricMatrix")
449		if(sym.x) {
450		W <- if(x@uplo == "U") # stored only [i,j] with i <= j
451		ij[,1] > ij[,2] else ij[,1] < ij[,2]
452		if(any(W))
453		ij[W,] <- ij[W, 2:1]
454		}
455		if(extends(cld, "sparseMatrix")) {
456		## do something smarter:
457		nr <- nrow(x)
458		if(!extends(cld, "CsparseMatrix")) {
459		x <- as(x, "CsparseMatrix") # simpler; our standard
460		cld <- getClassDef(class(x))
461		}
462		tri.x <- extends(cld, "triangularMatrix")
463		if(tri.x) {
464		## need these for the 'x' slot in any case
465		if (x@diag == "U") x <- .Call(Csparse_diagU2N, x)
466		## slightly more efficient than non0.i() or non0ind():
467		ij.x <- .Call(compressed_non_0_ij, x, isC=TRUE)
468		} else { ## symmetric / general : for symmetric, only "existing"b
469		ij.x <- non0.i(x, cld)
470		}
471
472		mi <- match(encodeInd(ij.x,   nr),
473		encodeInd(ij -1L, nr), nomatch=0)
474		mmi <- mi != 0
475		## Result:
476		ans <- vector(mode = .type.kind[.M.kindC(cld)], length = m)
477		## those that are *not* zero:
478		ans[mi[mmi]] <-
479		if(extends(cld, "nsparseMatrix")) TRUE else x@x[mmi]
480		ans
481
482		} else { ## non-sparse : dense
483		##---- NEVER happens:  'denseMatrix' has its own setMethod(.) !
484		message("m[ <ij-matrix> ]: inefficiently indexing single elements")
485		i1 <- ij[,1]
486		i2 <- ij[,2]
487		## very inefficient for large m
488		unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]]))
489		}
490		} else { # 1 <= m <= 3
491		i1 <- ij[,1]
492		i2 <- ij[,2]
493		unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]]))
494		}
495		}
496
497	## A[ ij ]  where ij is (i,j) 2-column matrix -- but also when that is logical mat!	## A[ ij ]  where ij is (i,j) 2-column matrix -- but also when that is logical mat!
498	.M.sub.i.2col <- function (x, i, j, ..., drop)	.M.sub.i.2col <- function (x, i, j, ..., drop)
499	{	{
#	Line 443	Line 509
509	m <- nrow(i)	m <- nrow(i)
510	if(m == 0) return(vector(mode = .type.kind[.M.kind(x)]))	if(m == 0) return(vector(mode = .type.kind[.M.kind(x)]))
511	## else	## else
512	i1 <- i[,1]	subset.ij(x, i)
	i2 <- i[,2]
	## potentially inefficient -- FIXME --
	unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]]))
513
514	} else stop("nargs() = ", nA,	} else stop("nargs() = ", nA,
515	".  Extraneous illegal arguments inside '[ .. ]' (i.2col)?")	".  Extraneous illegal arguments inside '[ .. ]' (i.2col)?")
#	Line 502	Line 565
565	value <- rep(value, length = m)	value <- rep(value, length = m)
566	i1 <- i[,1]	i1 <- i[,1]
567	i2 <- i[,2]	i2 <- i[,2]
568		if(m > 2)
569		message("m[ <ij-matrix> ] <- v: inefficiently treating single elements")
570	## inefficient -- FIXME -- (also loses "symmetry" unnecessarily)	## inefficient -- FIXME -- (also loses "symmetry" unnecessarily)
571	for(k in seq_len(m))	for(k in seq_len(m))
572	x[i1[k], i2[k]] <- value[k]	x[i1[k], i2[k]] <- value[k]

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.2110   changed lines   Added in v.2115

---

root@r-forge.r-project.org	ViewVC Help
Powered by ViewVC 1.0.0