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 1331, Sat Jul 22 17:59:53 2006 UTC revision 1617, Fri Oct 6 15:42:12 2006 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 16  Line 19 
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    ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
34    ### Bert's code built on a post by Andy Liaw who most probably was influenced
35    ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
36    ### who posted his bdiag() function written in December 1995.
37    
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  setAs("diagonalMatrix", "triangularMatrix",  setAs("diagonalMatrix", "triangularMatrix",
62        function(from) {        function(from) {
63            n <- from@Dim[1]            n <- from@Dim[1]
64            i <- seq(length = n)            i <- seq(length = n)
65            x <- from@x            x <- from@x
66            new(if(is.numeric(x)) "dtTMatrix" else "ltTMatrix",            new(paste(.M.kind(from), "tTMatrix", sep=''),
67                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
68                x = x, i = i, j = i)                x = x, i = i, j = i)
69            })            })
# Line 43  Line 77 
77        })        })
78    
79  setAs("diagonalMatrix", "generalMatrix",  setAs("diagonalMatrix", "generalMatrix",
80        function(from) {        function(from) as(from, paste(.M.kind(from), "geMatrix", sep='')))
           x <- as(from, "matrix")  
           as(x,  
              if(is.logical(x)) "lgeMatrix"  
 ## Not yet:  
 ##              else if(is.complex(x)) "zgeMatrix"  
 ##              else if(is.integer(x)) "igeMatrix"  
              else "dgeMatrix")  
       })  
81    
82  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
83        function(from) {        function(from) {
# Line 67  Line 93 
93  setAs("ldiMatrix", "lgTMatrix",  setAs("ldiMatrix", "lgTMatrix",
94        function(from) {        function(from) {
95            n <- from@Dim[1]            n <- from@Dim[1]
96            i <- (if(from@diag == "U") seq(length = n) else which(from@x)) - 1:1            if(from@diag == "U") { # unit-diagonal
97            new("lgTMatrix", i = i, j = i,                x <- rep.int(TRUE, n)
98                  i <- seq(length = n)
99              } else { # "normal"
100                  nz <- nz.NA(from@x, na. = TRUE)
101                  x <- from@x[nz]
102                  i <- which(nz) - 1:1
103              }
104              new("lgTMatrix", i = i, j = i, x = x,
105                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
106    
107  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
# Line 78  Line 111 
111        function(from)        function(from)
112            as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))            as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
113    
114    if(FALSE) # now have faster  "ddense" -> "dge"
115  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
116        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) as(as(from, "matrix"), "dgeMatrix"))
117    
# Line 95  Line 129 
129                cl <- "ddiMatrix"                cl <- "ddiMatrix"
130                uni <- all(x == 1)                uni <- all(x == 1)
131                storage.mode(x) <- "double"                storage.mode(x) <- "double"
132            }            } ## TODO: complex
133            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
134                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
135        })        })
# Line 120  Line 154 
154                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
155        })        })
156    
157    ## When you assign to a diagonalMatrix, the result should be
158    ## diagonal or sparse
159    setReplaceMethod("[", signature(x = "diagonalMatrix",
160                                    i = "ANY", j = "ANY", value = "ANY"),
161                     function(x, i, j, value) {
162                         r <- callGeneric(x = as(x,"sparseMatrix"),
163                                          i=i, j=j, value=value)
164                         if(isDiagonal(r)) as(r, "diagonalMatrix") else r
165                     })
166    
167    
168  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
169            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
170    
# Line 268  Line 313 
313  }  }
314    
315  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
316            function(object) prDiag(object))            function(object) {
317                  d <- dim(object)
318                  cl <- class(object)
319                  cat(sprintf('%d x %d diagonal matrix of class "%s"\n',
320                              d[1], d[2], cl))
321                  prDiag(object)
322              })

Legend:
Removed from v.1331  
changed lines
  Added in v.1617

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