SCM

SCM Repository

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

Diff of /pkg/Matrix/R/diagMatrix.R

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

revision 1575, Mon Sep 18 14:47:40 2006 UTC revision 1805, Tue Mar 27 16:46:03 2007 UTC
# Line 1  Line 1 
1    #### All methods for "diagonalMatrix" and its subclasses,
2    ####  currently "ddiMatrix", "ldiMatrix"
3    
4  ## Purpose: Constructor of diagonal matrices -- ~= diag() ,  ## Purpose: Constructor of diagonal matrices -- ~= diag() ,
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
# Line 10  Line 13 
13          n <- as.integer(n)          n <- as.integer(n)
14      }      }
15    
16      if(missing(x)) # unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          stopifnot(length(x) == n)          stopifnot(length(x) == n)
# Line 27  Line 30 
30      }      }
31  }  }
32    
33  setAs("diagonalMatrix", "triangularMatrix",  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
34        function(from) {  ### Bert's code built on a post by Andy Liaw who most probably was influenced
35            n <- from@Dim[1]  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
36            i <- seq(length = n)  ### who posted his bdiag() function written in December 1995.
37            x <- from@x  
38    bdiag <- function(...) {
39        if(nargs() == 0) return(new("dgCMatrix"))
40        ## else :
41        mlist <- if (nargs() == 1) as.list(...) else list(...)
42        dims <- sapply(mlist, dim)
43        ## make sure we had all matrices:
44        if(!(is.matrix(dims) && nrow(dims) == 2))
45            stop("some arguments are not matrices")
46        csdim <- rbind(rep.int(0:0, 2),
47                       apply(sapply(mlist, dim), 1, cumsum))
48        ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))
49        add1 <- matrix(1:0, 2,2)
50        for(i in seq_along(mlist)) {
51            indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
52            if(is.null(dim(indx))) ## non-square matrix
53                ret[indx[[1]],indx[[2]]] <- mlist[[i]]
54            else ## square matrix
55                ret[indx[,1],indx[,2]] <- mlist[[i]]
56        }
57        ## slightly debatable if we really should return Csparse.. :
58        as(ret, "CsparseMatrix")
59    }
60    
61    diag2T <- function(from) {
62        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1
63            new(paste(.M.kind(from), "tTMatrix", sep=''),            new(paste(.M.kind(from), "tTMatrix", sep=''),
64                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
65                x = x, i = i, j = i)          x = from@x, # <- ok for diag = "U" and "N" (!)
66            })          i = i, j = i)
67    }
68    
69    setAs("diagonalMatrix", "triangularMatrix", diag2T)
70    setAs("diagonalMatrix", "sparseMatrix", diag2T)
71    ## needed too (otherwise <dense> -> Tsparse is taken):
72    setAs("diagonalMatrix", "TsparseMatrix", diag2T)
73    ## is better than this:
74    ## setAs("diagonalMatrix", "sparseMatrix",
75    ##       function(from)
76    ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
77    setAs("diagonalMatrix", "CsparseMatrix",
78          function(from) as(diag2T(from), "CsparseMatrix"))
79    
80  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
81        function(from) {        function(from) {
# Line 45  Line 85 
85                 nrow = n, ncol = n)                 nrow = n, ncol = n)
86        })        })
87    
88  setAs("diagonalMatrix", "generalMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
89        function(from) as(from, paste(.M.kind(from), "geMatrix", sep='')))        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))
90    
91    .diag.x <- function(m) {
92        if(m@diag == "U")
93            rep.int(if(is.numeric(m@x)) 1. else TRUE,
94                    m@Dim[1])
95        else m@x
96    }
97    
98    .diag.2N <- function(m) {
99        if(m@diag == "U") m@diag <- "N"
100        m
101    }
102    
103    ## given the above, the following  4  coercions should be all unneeded;
104    ## we prefer triangular to general:
105  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
106        function(from) {        function(from) {
107              .Deprecated("as(, \"sparseMatrix\")")
108            n <- from@Dim[1]            n <- from@Dim[1]
109            i <- seq(length = n) - 1:1            i <- seq_len(n) - 1:1
110            new("dgTMatrix", i = i, j = i,            new("dgTMatrix", i = i, j = i, x = .diag.x(from),
               x = if(from@diag == "U") rep(1,n) else from@x,  
111                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
112    
113  setAs("ddiMatrix", "dgCMatrix",  setAs("ddiMatrix", "dgCMatrix",
114        function(from) as(as(from, "dgTMatrix"), "dgCMatrix"))        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))
115    
116  setAs("ldiMatrix", "lgTMatrix",  setAs("ldiMatrix", "lgTMatrix",
117        function(from) {        function(from) {
118              .Deprecated("as(, \"sparseMatrix\")")
119            n <- from@Dim[1]            n <- from@Dim[1]
120            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
121                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
122                i <- seq(length = n)                i <- seq_len(n) - 1:1
123            } else { # "normal"            } else { # "normal"
124                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
125                x <- from@x[nz]                x <- from@x[nz]
# Line 76  Line 131 
131  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
132        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
133    
 setAs("diagonalMatrix", "sparseMatrix",  
       function(from)  
           as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  
134    
135  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
136  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
# Line 123  Line 175 
175                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
176        })        })
177    
178    
179    setMethod("diag", signature(x = "diagonalMatrix"),
180              function(x = 1, nrow, ncol = n) .diag.x(x))
181    
182    
183    subDiag <- function(x, i, j, drop) {
184        x <- as(x, "sparseMatrix")
185        x <- if(missing(i))
186            x[, j, drop=drop]
187        else if(missing(j))
188            x[i, , drop=drop]
189        else
190            x[i,j, drop=drop]
191        if(isDiagonal(x)) as(x, "diagonalMatrix") else x
192    }
193    
194    setMethod("[", signature(x = "diagonalMatrix", i = "index",
195                             j = "index", drop = "logical"), subDiag)
196    setMethod("[", signature(x = "diagonalMatrix", i = "index",
197                            j = "missing", drop = "logical"),
198              function(x, i, drop) subDiag(x, i=i, drop=drop))
199    setMethod("[", signature(x = "diagonalMatrix", i = "missing",
200                             j = "index", drop = "logical"),
201              function(x, j, drop) subDiag(x, j=j, drop=drop))
202    
203    ## When you assign to a diagonalMatrix, the result should be
204    ## diagonal or sparse ---
205    ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
206    replDiag <- function(x, i, j, value) {
207        x <- as(x, "sparseMatrix")
208        if(missing(i))
209            x[, j] <- value
210        else if(missing(j))
211            x[i, ] <- value
212        else
213            x[i,j] <- value
214        if(isDiagonal(x)) as(x, "diagonalMatrix") else x
215    }
216    
217    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
218                                    j = "index", value = "replValue"), replDiag)
219    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
220                                    j = "missing", value = "replValue"),
221                     function(x, i, value) replDiag(x, i=i, value=value))
222    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
223                                    j = "index", value = "replValue"),
224                     function(x, j, value) replDiag(x, j=j, value=value))
225    
226    
227  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
228            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
229    
# Line 133  Line 234 
234  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
235            function(object) TRUE)            function(object) TRUE)
236    
237  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
238            function(x = 1, nrow, ncol = n) {            function(x, pivot) {
239               if(x@diag == "U")                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")
240                   rep.int(if(is.logical(x@x)) TRUE else 1, x@Dim[1])                x@x <- sqrt(x@x)
241               else x@x                x
242            })            })
243    ## chol(L) is L for logical diagonal:
244    setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)
245    
246  setMethod("!", "ldiMatrix", function(e1) {  setMethod("!", "ldiMatrix", function(e1) {
247      if(e1@diag == "N")      if(e1@diag == "N")
# Line 147  Line 250 
250          e1@diag <- "N"          e1@diag <- "N"
251          e1@x <- rep.int(FALSE, e1@Dim[1])          e1@x <- rep.int(FALSE, e1@Dim[1])
252      }      }
253      x      e1
254  })  })
255    
256  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
257    ##       ---------------------
258  ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"  ## Note that "ldi" logical are treated as numeric
259  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
260      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      if(any(dim(x) != dim(y))) stop("non-matching dimensions")
261      if(x@diag != "U") {      if(x@diag != "U") {
262          if(y@diag != "U") x@x <- x@x * y@x          if(y@diag != "U") {
263                nx <- x@x * y@x
264                if(is.numeric(nx) && !is.numeric(x@x))
265                    x <- as(x, "dMatrix")
266                x@x <- as.numeric(nx)
267            }
268          return(x)          return(x)
269      } else ## x is unit diagonal      } else ## x is unit diagonal
270      return(y)      return(y)
271  }  }
272    
273  setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
274            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
275    
276  formals(diagdiagprod) <- alist(x=, y=NULL)  formals(diagdiagprod) <- alist(x=, y=x)
277  setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
278              diagdiagprod, valueClass = "ddiMatrix")
279    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
280              diagdiagprod, valueClass = "ddiMatrix")
281    setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
282            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
283  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
284            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
285    
286    
# Line 185  Line 297 
297  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
298  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
299            diagmatprod)            diagmatprod)
300    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),
301              diagmatprod)
302    
303  diagdgeprod <- function(x, y) {  diagdgeprod <- function(x, y) {
304      dx <- dim(x)      dx <- dim(x)
# Line 205  Line 319 
319                dx <- dim(x)                dx <- dim(x)
320                dy <- dim(y)                dy <- dim(y)
321                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
322                as(if(y@diag == "U") x else x * rep.int(y@x, dx[1]), "Matrix")                as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix")
323            })            })
324    
325  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),
# Line 214  Line 328 
328                dy <- dim(y)                dy <- dim(y)
329                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
330                if(y@diag == "N")                if(y@diag == "N")
331                    x@x <- x@x * rep.int(y@x, dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
332                x                x
333            })            })
334    
# Line 231  Line 345 
345  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
346  ##           })  ##           })
347    
348    setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
349              function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })
350    
351  ### ---------------- diagonal  o   sparse  -----------------------------  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
352              function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })
353    
354  ## These are cheap implementations via coercion  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
355              function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })
356    
357    setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
358              function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })
359    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
360    
361    ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
362  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
363            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "sparseMatrix") %*% y)
364    ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
365    ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
366    ## ==> do this:
367    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
368              function(x, y) as(x, "CsparseMatrix") %*% y)
369    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
370              function(x, y) x %*% as(y, "CsparseMatrix"))
371    ## NB: this is *not* needed for Tsparse & Rsparse
372    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
373    ##       do indeed work by going through sparse (and *not* ddense)!
374    
375  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
376            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) x %*% as(y, "sparseMatrix"))
377    
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  
           function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })  
378    
379  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
380            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(a, b, ...) {
381                  a@x <- 1/ a@x
382                  a@Dimnames <- a@Dimnames[2:1]
383                  a
384              })
385    
386    solveDiag <- function(a, b, ...) {
387        if((n <- a@Dim[1]) != nrow(b))
388            stop("incompatible matrix dimensions")
389        ## trivially invert a 'in place' and multiply:
390        a@x <- 1/ a@x
391        a@Dimnames <- a@Dimnames[2:1]
392        a %*% b
393    }
394    setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"),
395              solveDiag)
396    setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
397              solveDiag)
398    
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  
           function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })  
399    
 setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })  
400    
401    
402    ### ---------------- diagonal  o  sparse  -----------------------------
403    
404    
405    ## Use function for several signatures, in order to evade
406    ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
407    diagOdiag <- function(e1,e2) { # result should also be diagonal
408        r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
409        if(is.numeric(r)) {
410            if(is.numeric(e2@x)) {
411                e2@x <- r; return(.diag.2N(e2)) }
412            if(!is.numeric(e1@x))
413                ## e.g. e1, e2 are logical;
414                e1 <- as(e1, "dMatrix")
415        }
416        else if(is.logical(r))
417            e1 <- as(e1, "lMatrix")
418        else stop("intermediate 'r' is of type", typeof(r))
419        e1@x <- r
420        .diag.2N(e1)
421    }
422    
423    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
424              diagOdiag)
425    ## These two are just for method disambiguation:
426    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),
427              diagOdiag)
428    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
429              diagOdiag)
430    
431    ## For almost everything else, diag* shall be treated "as sparse" :
432    ## These are cheap implementations via coercion
433    
434    ## for disambiguation
435    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
436              function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))
437    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
438              function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))
439    ## in general:
440    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
441              function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))
442    setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
443              function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))
444    
445    
446    
447  ## similar to prTriang() in ./Auxiliaries.R :  ## similar to prTriang() in ./Auxiliaries.R :
# Line 271  Line 456 
456  }  }
457    
458  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
459            function(object) prDiag(object))            function(object) {
460                  d <- dim(object)
461                  cl <- class(object)
462                  cat(sprintf('%d x %d diagonal matrix of class "%s"\n',
463                              d[1], d[2], cl))
464                  prDiag(object)
465              })

Legend:
Removed from v.1575  
changed lines
  Added in v.1805

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