SCM

SCM Repository

[matrix] Diff of /pkg/R/Matrix.R
ViewVC logotype

Diff of /pkg/R/Matrix.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 871, Fri Aug 26 17:26:49 2005 UTC revision 1315, Mon Jul 10 16:08:05 2006 UTC
# Line 1  Line 1 
1  #### Toplevel ``virtual'' class "Matrix"  #### Toplevel ``virtual'' class "Matrix"
2    
3  ## probably not needed eventually:  ## ## probably not needed eventually:
4  setAs(from = "ddenseMatrix", to = "matrix",  ## setAs(from = "ddenseMatrix", to = "matrix",
5        function(from) {  ##       function(from) {
6            if(length(d <- dim(from)) != 2) stop("dim(.) has not length 2")  ##        if(length(d <- dim(from)) != 2) stop("dim(.) has not length 2")
7            array(from@x, dim = d, dimnames = dimnames(from))  ##        array(from@x, dim = d, dimnames = dimnames(from))
8        })  ##       })
   
 ## private function to be used as show() method possibly more than once  
 prMatrix <- function(object) {  
     d <- dim(object)  
     cl <- class(object)  
     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))  
     m <- as(object, "matrix")  
     maxp <- getOption("max.print")  
     if(prod(d) <= maxp) print(m)  
     else { ## d[1] > maxp / d[2] >= nr :  
         nr <- maxp %/% d[2]  
         n2 <- ceiling(nr / 2)  
         print(head(m, max(1, n2)))  
         cat("\n ..........\n\n")  
         print(tail(m, max(1, nr - n2)))  
     }  
     ## DEBUG: cat("str(.):\n") ; str(object)  
     invisible(object)# as print() S3 methods do  
 }  
   
 setMethod("show", signature(object = "ddenseMatrix"), prMatrix)  
   
 ##- ## FIXME: The following is only for the "dMatrix" objects that are not  
 ##- ##        "dense" nor "sparse" -- i.e. "packed" ones :  
 ##- ## But these could be printed better -- "." for structural zeros.  
 ##- setMethod("show", signature(object = "dMatrix"), prMatrix)  
 ##- ## and improve this as well:  
 ##- setMethod("show", signature(object = "pMatrix"), prMatrix)  
 ##- ## this should now be superfluous [keep for safety for the moment]:  
 setMethod("show", signature(object = "Matrix"), prMatrix)  
9    
10  ## should propagate to all subclasses:  ## should propagate to all subclasses:
11  setMethod("as.matrix", signature(x = "Matrix"), function(x) as(x, "matrix"))  setMethod("as.matrix", signature(x = "Matrix"), function(x) as(x, "matrix"))
12    ## for 'Matrix' objects, as.array() should be equivalent:
13    setMethod("as.array",  signature(x = "Matrix"), function(x) as(x, "matrix"))
14    
15    ## slow "fall back" method {subclasses should have faster ones}:
16    setMethod("as.vector", signature(x = "Matrix", mode = "missing"),
17              function(x) as.vector(as(x, "matrix")))
18    
19    
20    ## Note that isSymmetric is *not* exported
21    ## but that "base" has an isSymmetric() S3-generic since R 2.3.0
22    setMethod("isSymmetric", signature(object = "symmetricMatrix"),
23              function(object,tol) TRUE)
24    setMethod("isSymmetric", signature(object = "triangularMatrix"),
25              ## TRUE iff diagonal:
26              function(object,tol) isDiagonal(object))
27    
28    if(paste(R.version$major, R.version$minor, sep=".") < "2.3")
29        ## need a "matrix" method as in R 2.3 and later
30        setMethod("isSymmetric", signature(object = "matrix"),
31                  function(object, tol = 100*.Machine$double.eps, ...)
32              {
33                  ## pretest: is it square?
34                  d <- dim(object)
35                  if(d[1] != d[2]) return(FALSE)
36                  ## for `broken' all.equal in R <= 2.2.x:
37                  dn <- dimnames(object)
38                  if(!identical(dn[1], dn[2])) return(FALSE)
39                  test <-
40                      if(is.complex(object))
41                          all.equal.numeric(object, Conj(t(object)), tol = tol, ...)
42                      else              # numeric, character, ..
43                          all.equal(object, t(object), tol = tol, ...)
44                  isTRUE(test)
45              })
46    
47    
48    setMethod("isTriangular", signature(object = "triangularMatrix"),
49              function(object, ...) TRUE)
50    
51    setMethod("isTriangular", signature(object = "matrix"), isTriMat)
52    
53    setMethod("isDiagonal", signature(object = "matrix"), .is.diagonal)
54    
55    
56    
57  setMethod("dim", signature(x = "Matrix"),  setMethod("dim", signature(x = "Matrix"),
58            function(x) x@Dim, valueClass = "integer")            function(x) x@Dim, valueClass = "integer")
# Line 62  Line 76 
76            function(obj) { obj@Dimnames <- list(NULL,NULL); obj})            function(obj) { obj@Dimnames <- list(NULL,NULL); obj})
77    
78  Matrix <-  Matrix <-
79      function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)      function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL,
80                  sparse = NULL)
81  {  {
82      if (is(data, "Matrix")) return(data)      sparseDefault <- function(m)
83      if (is.matrix(data)) { val <- data }          prod(dim(m)) > 2*sum(as(m, "matrix") != 0)
84      else { ## cut & paste from "base::matrix" :  
85        i.M <- is(data, "Matrix")
86        if(is.null(sparse) && (i.M || is(data, "matrix")))
87            sparse <- sparseDefault(data)
88    
89        if (i.M) {
90            sM <- is(data,"sparseMatrix")
91            if((sparse && sM) || (!sparse && !sM))
92                return(data)
93            ## else : convert  dense <-> sparse -> at end
94        }
95        else if (!is.matrix(data)) { ## cut & paste from "base::matrix" :
96          if (missing(nrow))          if (missing(nrow))
97              nrow <- ceiling(length(data)/ncol)              nrow <- ceiling(length(data)/ncol)
98          else if (missing(ncol))          else if (missing(ncol))
99              ncol <- ceiling(length(data)/nrow)              ncol <- ceiling(length(data)/nrow)
100          val <- .Internal(matrix(data, nrow, ncol, byrow))          if(length(data) == 1 && data == 0 && !identical(sparse,FALSE)) {
101          dimnames(val) <- dimnames              if(is.null(sparse)) sparse <- TRUE
102                ## will be sparse: do NOT construct full matrix!
103                data <- new(if(is.numeric(data)) "dgTMatrix" else
104                            if(is.logical(data)) "lgTMatrix" else
105                            stop("invalid 'data'"),
106                            Dim = as.integer(c(nrow,ncol)),
107                            Dimnames = if(is.null(dimnames))
108                            list(NULL,NULL) else dimnames)
109            } else { ## normal case
110                data <- .Internal(matrix(data, nrow, ncol, byrow))
111                if(is.null(sparse))
112                    sparse <- sparseDefault(data)
113                dimnames(data) <- dimnames
114            }
115        } else if (!is.null(dimnames))
116            dimnames(data) <- dimnames
117    
118        ## 'data' is now a "matrix" or "Matrix"
119    
120        ## check for symmetric / triangular / diagonal :
121        isSym <- isSymmetric(data)
122        if((isTri <- !isSym))
123            isTri <- isTriangular(data)
124        isDiag <- isSym # cannot be diagonal if it isn't symmetric
125        if(isDiag)
126            isDiag <- isDiagonal(data)
127    
128    ### TODO: Compare with as.Matrix() and its tests in ./dgeMatrix.R
129    
130        ## Find proper matrix class 'cl'
131        cl <-
132            if(isDiag)
133                "diagonalMatrix" # -> will automatically check for type
134            else {
135                ## consider it's type
136                ctype <-
137                    if(is(data,"Matrix")) class(data)
138                    else {
139                        if("complex" == (ctype <- typeof(data)))
140                            "z" else ctype
141                    }
142                ctype <- substr(ctype, 1,1) # "d", "l", "i" or "z"
143                if(ctype == "z")
144                    stop("complex matrices not yet implemented in Matrix package")
145                if(ctype == "i") {
146                    warning("integer matrices not yet implemented in 'Matrix'; ",
147                            "using 'double' ones'")
148                    ctype <- "d"
149      }      }
150      as(val, "dgeMatrix")              paste(ctype,
151                      if(sparse) {
152                          if(isSym) "sCMatrix" else
153                          if(isTri) "tCMatrix" else "gCMatrix"
154                      } else { ## dense
155                          if(isSym) "syMatrix" else
156                          if(isTri) "trMatrix" else "geMatrix"
157                      }, sep="")
158            }
159    
160        ## Now coerce and return
161        as(data, cl)
162  }  }
163    
164  ## Methods for operations where one argument is numeric  ## Methods for operations where one argument is numeric
# Line 97  Line 181 
181  setMethod("solve", signature(a = "Matrix", b = "numeric"),  setMethod("solve", signature(a = "Matrix", b = "numeric"),
182            function(a, b, ...) callGeneric(a, as.matrix(b)))            function(a, b, ...) callGeneric(a, as.matrix(b)))
183    
184    ## bail-out methods in order to get better error messages
185    setMethod("%*%", signature(x = "Matrix", y = "Matrix"),
186              function (x, y)
187              stop(gettextf('not-yet-implemented method for <%s> %%*%% <%s>',
188                            class(x), class(y))))
189    
190    setMethod("crossprod", signature(x = "Matrix", y = "ANY"),
191              function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y)))
192    setMethod("crossprod", signature(x = "ANY", y = "Matrix"),
193              function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y)))
194    
195    ## There are special sparse methods; this is a "fall back":
196    setMethod("kronecker", signature(X = "Matrix", Y = "ANY",
197                                     FUN = "ANY", make.dimnames = "ANY"),
198              function(X, Y, FUN, make.dimnames, ...) {
199                  X <- as(X, "matrix") ; Matrix(callGeneric()) })
200    setMethod("kronecker", signature(X = "ANY", Y = "Matrix",
201                                     FUN = "ANY", make.dimnames = "ANY"),
202              function(X, Y, FUN, make.dimnames, ...) {
203                  Y <- as(Y, "matrix") ; Matrix(callGeneric()) })
204    
205    
206    setMethod("diag", signature(x = "Matrix"),
207              function(x, nrow, ncol) .bail.out.1(.Generic, class(x)))
208    setMethod("t", signature(x = "Matrix"),
209              function(x) .bail.out.1(.Generic, class(x)))
210    
211    ## Group Methods
212    setMethod("+", signature(e1 = "Matrix", e2 = "missing"), function(e1) e1)
213    ## "fallback":
214    setMethod("-", signature(e1 = "Matrix", e2 = "missing"),
215              function(e1) {
216                  warning("inefficient method used for \"- e1\"")
217                  0-e1
218              })
219    
220    ## bail-outs:
221    setMethod("Compare", signature(e1 = "Matrix", e2 = "Matrix"),
222              function(e1, e2) {
223                  d <- dimCheck(e1,e2)
224                  .bail.out.2(.Generic, class(e1), class(e2))
225              })
226    setMethod("Compare", signature(e1 = "Matrix", e2 = "ANY"),
227              function(e1, e2) .bail.out.2(.Generic, class(e1), class(e2)))
228    setMethod("Compare", signature(e1 = "ANY", e2 = "Matrix"),
229              function(e1, e2) .bail.out.2(.Generic, class(e1), class(e2)))
230    
231    
232    
233  ### --------------------------------------------------------------------------  ### --------------------------------------------------------------------------
234  ###  ###
235  ### Subsetting "["  and  ### Subsetting "["  and
236  ### SubAssign  "[<-" : The "missing" cases can be dealt with here, "at the top":  ### SubAssign  "[<-" : The "missing" cases can be dealt with here, "at the top":
237    
238    ## Using "index" for indices should allow
239    ## integer (numeric), logical, or character (names!) indices :
240    
241  ## "x[]":  ## "x[]":
242  setMethod("[", signature(x = "Matrix",  setMethod("[", signature(x = "Matrix",
243                           i = "missing", j = "missing", drop = "ANY"),                           i = "missing", j = "missing", drop = "ANY"),
244            function (x, i, j, drop) x)            function (x, i, j, drop) x)
245    
246  ## missing 'drop' --> 'drop = TRUE'  ## missing 'drop' --> 'drop = TRUE'
247  ##                     -----------  ##                     -----------
248  ## select rows  ## select rows
249  setMethod("[", signature(x = "Matrix", i = "numeric", j = "missing",  setMethod("[", signature(x = "Matrix", i = "index", j = "missing",
250                           drop = "missing"),                           drop = "missing"),
251            function(x,i,j, drop) callGeneric(x, i=i, drop= TRUE))            function(x,i,j, drop) callGeneric(x, i=i, drop= TRUE))
252  ## select columns  ## select columns
253  setMethod("[", signature(x = "Matrix", i = "missing", j = "numeric",  setMethod("[", signature(x = "Matrix", i = "missing", j = "index",
254                           drop = "missing"),                           drop = "missing"),
255            function(x,i,j, drop) callGeneric(x, j=j, drop= TRUE))            function(x,i,j, drop) callGeneric(x, j=j, drop= TRUE))
256  setMethod("[", signature(x = "Matrix", i = "numeric", j = "numeric",  setMethod("[", signature(x = "Matrix", i = "index", j = "index",
257                           drop = "missing"),                           drop = "missing"),
258            function(x,i,j, drop) callGeneric(x, i=i, j=j, drop= TRUE))            function(x,i,j, drop) callGeneric(x, i=i, j=j, drop= TRUE))
259    
260    ## bail out if any of (i,j,drop) is "non-sense"
261    setMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", drop = "ANY"),
262              function(x,i,j, drop)
263              stop("invalid or not-yet-implemented 'Matrix' subsetting"))
264    
265    ##  "logical *vector* indexing, such as  M [ M >= 10 ] :
266    setMethod("[", signature(x = "Matrix", i = "lMatrix", j = "missing",
267                             drop = "ANY"),
268              function (x, i, j, drop) {
269                  as(x, geClass(x))@x[as.vector(i)]
270                                            # -> error when lengths don't match
271              })
272    
273    setMethod("[", signature(x = "Matrix", i = "logical", j = "missing",
274                             drop = "ANY"),
275              function (x, i, j, drop) as(x, geClass(x))@x[i])
276    
277    
278  ## "FIXME:"  ## "FIXME:"
279  ## How can we get at   A[ ij ]  where ij is (i,j) 2-column matrix?  ## How can we get at   A[ ij ]  where ij is (i,j) 2-column matrix?
280  ##  and                A[ LL ]  where LL is a logical *vector*  ##  and                A[ LL ]  where LL is a logical *vector*
281    ## -> [.data.frame uses nargs() - can we do this in the *generic* ?
282    
283    
284  ### "[<-" : -----------------  ### "[<-" : -----------------
285    
286  ## x[] <- value :  ## x[] <- value :
287  setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "missing",  setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "missing",
288                                  value = "vector"),##  double/logical/...                                  value = "ANY"),## double/logical/...
289            function (x, value) { x@x <- value ; validObject(x); x })            function (x, value) {
290                  x@x <- value
291                  validObject(x)# check if type and lengths above match
292                  x
293              })
294    
295  ## Otherwise (value is not "vector"): bail out  ## Method for all 'Matrix' kinds (rather than incomprehensible error messages);
296    ## (ANY,ANY,ANY) is used when no `real method' is implemented :
297  setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY",  setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY",
298                                  value = "ANY"),                                  value = "ANY"),
299            function (x, i, j, value) stop("RHS 'value' must be of class \"vector\""))            function (x, i, j, value) {
300                  if(!is.atomic(value))
301                      stop("RHS 'value' must match matrix class ", class(x))
302                  else stop("not-yet-implemented 'Matrix[<-' method")
303              })
304    
305    
306    ## The trivial methods :
307    setMethod("cbind2", signature(x = "Matrix", y = "NULL"),
308              function(x, y) x)
309    setMethod("cbind2", signature(x = "Matrix", y = "missing"),
310              function(x, y) x)
311    setMethod("cbind2", signature(x = "NULL", y="Matrix"),
312              function(x, y) x)
313    
314    setMethod("rbind2", signature(x = "Matrix", y = "NULL"),
315              function(x, y) x)
316    setMethod("rbind2", signature(x = "Matrix", y = "missing"),
317              function(x, y) x)
318    setMethod("rbind2", signature(x = "NULL", y="Matrix"),
319              function(x, y) x)
320    
321    ## Makes sure one gets x decent error message for the unimplemented cases:
322    setMethod("cbind2", signature(x = "Matrix", y = "Matrix"),
323              function(x, y) {
324                  rowCheck(x,y)
325                  stop(gettextf("cbind2() method for (%s,%s) not-yet defined",
326                                class(x), class(y)))
327              })
328    
329  if(FALSE) ## The following can't work as long as cbind is function(..., *)  ## Use a working fall back {particularly useful for sparse}:
330  setMethod("cbind", signature(a = "Matrix", b = "Matrix"),  ## FIXME: implement rbind2 via "cholmod" for C* and Tsparse ones
331            function(a, b, ...) {  setMethod("rbind2", signature(x = "Matrix", y = "Matrix"),
332                da <- Dim(a)            function(x, y) {
333                db <- Dim(b)                colCheck(x,y)
334                if(da[1] != db[1])                t(cbind2(t(x), t(y)))
                   stop("Matrices must have same number of rows for cbind()ing")  
335            })            })

Legend:
Removed from v.871  
changed lines
  Added in v.1315

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge