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 1617, Fri Oct 6 15:42:12 2006 UTC revision 1805, Tue Mar 27 16:46:03 2007 UTC
# Line 13  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 47  Line 47 
47                     apply(sapply(mlist, dim), 1, cumsum))                     apply(sapply(mlist, dim), 1, cumsum))
48      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))
49      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
50      for(i in seq(along = mlist)) {      for(i in seq_along(mlist)) {
51          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
52          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
53              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              ret[indx[[1]],indx[[2]]] <- mlist[[i]]
# Line 58  Line 58 
58      as(ret, "CsparseMatrix")      as(ret, "CsparseMatrix")
59  }  }
60    
61  setAs("diagonalMatrix", "triangularMatrix",  diag2T <- function(from) {
62        function(from) {      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1
           n <- from@Dim[1]  
           i <- seq(length = n)  
           x <- from@x  
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 76  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 107  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 154  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  ## When you assign to a diagonalMatrix, the result should be
204  ## diagonal or sparse  ## diagonal or sparse ---
205  setReplaceMethod("[", signature(x = "diagonalMatrix",  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
206                                  i = "ANY", j = "ANY", value = "ANY"),  replDiag <- function(x, i, j, value) {
207                   function(x, i, j, value) {      x <- as(x, "sparseMatrix")
208                       r <- callGeneric(x = as(x,"sparseMatrix"),      if(missing(i))
209                                        i=i, j=j, value=value)          x[, j] <- value
210                       if(isDiagonal(r)) as(r, "diagonalMatrix") else r      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"),
# Line 175  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 189  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")            diagdiagprod, valueClass = "ddiMatrix")
281  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  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 227  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 247  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 256  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 273  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  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  solveDiag <- function(a, b, ...) {
387            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })      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 = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })  
399    
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 :
448  prDiag <-  prDiag <-

Legend:
Removed from v.1617  
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