SCM

SCM Repository

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

Diff of /pkg/R/diagMatrix.R

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

revision 1635, Tue Oct 17 21:21:44 2006 UTC revision 1655, Mon Oct 30 17:16:27 2006 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    ## is better than this:
72    ## setAs("diagonalMatrix", "sparseMatrix",
73    ##       function(from)
74    ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
75    setAs("diagonalMatrix", "CsparseMatrix",
76          function(from) as(diag2T(from), "CsparseMatrix"))
77    
78  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
79        function(from) {        function(from) {
# Line 76  Line 83 
83                 nrow = n, ncol = n)                 nrow = n, ncol = n)
84        })        })
85    
86  setAs("diagonalMatrix", "generalMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
87        function(from) as(from, paste(.M.kind(from), "geMatrix", sep='')))        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))
88    
89    .diag.x <- function(m) {
90        if(m@diag == "U")
91            rep.int(if(is.numeric(m@x)) 1. else TRUE,
92                    m@Dim[1])
93        else m@x
94    }
95    
96    .diag.2N <- function(m) {
97        if(m@diag == "U") m@diag <- "N"
98        m
99    }
100    
101    ## given the above, the following  4  coercions should be all unneeded;
102    ## we prefer triangular to general:
103  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
104        function(from) {        function(from) {
105              .Deprecated("as(, \"sparseMatrix\")")
106            n <- from@Dim[1]            n <- from@Dim[1]
107            i <- seq(length = n) - 1:1            i <- seq_len(n) - 1:1
108            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,  
109                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
110    
111  setAs("ddiMatrix", "dgCMatrix",  setAs("ddiMatrix", "dgCMatrix",
112        function(from) as(as(from, "dgTMatrix"), "dgCMatrix"))        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))
113    
114  setAs("ldiMatrix", "lgTMatrix",  setAs("ldiMatrix", "lgTMatrix",
115        function(from) {        function(from) {
116              .Deprecated("as(, \"sparseMatrix\")")
117            n <- from@Dim[1]            n <- from@Dim[1]
118            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
119                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
120                i <- seq(length = n)                i <- seq_len(n) - 1:1
121            } else { # "normal"            } else { # "normal"
122                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
123                x <- from@x[nz]                x <- from@x[nz]
# Line 107  Line 129 
129  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
130        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
131    
 setAs("diagonalMatrix", "sparseMatrix",  
       function(from)  
           as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  
132    
133  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
134  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
# Line 175  Line 194 
194  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
195            function(object) TRUE)            function(object) TRUE)
196    
197    setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
198              function(x, pivot) {
199                  if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")
200                  x@x <- sqrt(x@x)
201                  x
202              })
203    ## chol(L) is L for logical diagonal:
204    setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)
205    
206    
207  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
208            function(x = 1, nrow, ncol = n) {            function(x = 1, nrow, ncol = n) {
209               if(x@diag == "U")               if(x@diag == "U")
# Line 193  Line 222 
222  })  })
223    
224  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
225    ##       ---------------------
226  ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"  ## Note that "ldi" logical are treated as numeric
227  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
228      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      if(any(dim(x) != dim(y))) stop("non-matching dimensions")
229      if(x@diag != "U") {      if(x@diag != "U") {
230          if(y@diag != "U") x@x <- x@x * y@x          if(y@diag != "U") {
231                nx <- x@x * y@x
232                if(is.numeric(nx) && !is.numeric(x@x))
233                    x <- as(x, "dMatrix")
234                x@x <- as.numeric(nx)
235            }
236          return(x)          return(x)
237      } else ## x is unit diagonal      } else ## x is unit diagonal
238      return(y)      return(y)
239  }  }
240    
241  setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
242            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
243    
244  formals(diagdiagprod) <- alist(x=, y=NULL)  formals(diagdiagprod) <- alist(x=, y=x)
245  setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
246            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
247  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
248              diagdiagprod, valueClass = "ddiMatrix")
249    setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
250              diagdiagprod, valueClass = "ddiMatrix")
251    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
252            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
253    
254    
# Line 227  Line 265 
265  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
266  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
267            diagmatprod)            diagmatprod)
268    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),
269              diagmatprod)
270    
271  diagdgeprod <- function(x, y) {  diagdgeprod <- function(x, y) {
272      dx <- dim(x)      dx <- dim(x)
# Line 276  Line 316 
316    
317  ### ---------------- diagonal  o   sparse  -----------------------------  ### ---------------- diagonal  o   sparse  -----------------------------
318    
319    
320    ## Use function for several signatures, in order to evade
321    ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
322    diagOdiag <- function(e1,e2) { # result should also be diagonal
323        r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
324        if(is.numeric(r)) {
325            if(is.numeric(e2@x)) {
326                e2@x <- r; return(.diag.2N(e2)) }
327            if(!is.numeric(e1@x))
328                ## e.g. e1, e2 are logical;
329                e1 <- as(e1, "dMatrix")
330        }
331        else if(is.logical(r))
332            e1 <- as(e1, "lMatrix")
333        else stop("intermediate 'r' is of type", typeof(r))
334        e1@x <- r
335        .diag.2N(e1)
336    }
337    
338    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
339              diagOdiag)
340    ## These two are just for method disambiguation:
341    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),
342              diagOdiag)
343    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
344              diagOdiag)
345    
346    ## For almost everything else, diag* shall be treated "as sparse" :
347  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
348    
349    ## for disambiguation
350    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
351              function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))
352    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
353              function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))
354    ## in general:
355    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
356              function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))
357    setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
358              function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))
359    
360    
361    
362  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
363    
364  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),

Legend:
Removed from v.1635  
changed lines
  Added in v.1655

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