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 1447, Fri Aug 25 20:38:36 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)
20          if(is.logical(x))          if(is.logical(x))
21              cl <- "ldiMatrix"              cl <- "ldiMatrix"
22          else {          else if(is.numeric(x)) {
23              cl <- "ddiMatrix"              cl <- "ddiMatrix"
24              x <- as.numeric(x)              x <- as.numeric(x)
25          }          }
26            else if(is.complex(x)) {
27                cl <- "zdiMatrix"  # will not yet work
28            } else stop("'x' has invalid data type")
29          new(cl, Dim = c(n,n), diag = "N", x = x)          new(cl, Dim = c(n,n), diag = "N", x = x)
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            new(if(is.numeric(x)) "dtTMatrix" else "ltTMatrix",  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=''),
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 42  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) {        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))
90            x <- as(from, "matrix")  
91            as(x,  .diag.x <- function(m) {
92               if(is.logical(x)) "lgeMatrix"      if(m@diag == "U")
93  ## Not yet:          rep.int(if(is.numeric(m@x)) 1. else TRUE,
94  ##              else if(is.complex(x)) "zgeMatrix"                  m@Dim[1])
95  ##              else if(is.integer(x)) "igeMatrix"      else m@x
96               else "dgeMatrix")  }
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            i <- (if(from@diag == "U") seq(length = n) else which(from@x)) - 1:1            if(from@diag == "U") { # unit-diagonal
121            new("lgTMatrix", i = i, j = i,                x <- rep.int(TRUE, n)
122                  i <- seq_len(n) - 1:1
123              } else { # "normal"
124                  nz <- nz.NA(from@x, na. = TRUE)
125                  x <- from@x[nz]
126                  i <- which(nz) - 1:1
127              }
128              new("lgTMatrix", i = i, j = i, x = x,
129                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
130    
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 96  Line 150 
150                cl <- "ddiMatrix"                cl <- "ddiMatrix"
151                uni <- all(x == 1)                uni <- all(x == 1)
152                storage.mode(x) <- "double"                storage.mode(x) <- "double"
153            }            } ## TODO: complex
154            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
155                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
156        })        })
# Line 121  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 131  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 145  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")            diagdiagprod, valueClass = "ddiMatrix")
279  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
280              diagdiagprod, valueClass = "ddiMatrix")
281    setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
282              diagdiagprod, valueClass = "ddiMatrix")
283    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
284            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
285    
286    
# Line 183  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 203  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 212  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 229  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  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
358              function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })
359    
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 269  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.1447  
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