# SCM Repository

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

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

revision 2814, Tue Jul 24 14:02:28 2012 UTC revision 2820, Mon Aug 20 14:06:23 2012 UTC
# Line 5  Line 5
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      if(missing(n))
10          n <- length(x)          n <- length(x)
11      else {      else {
# Line 17  Line 17
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          lx <- length(x)          lx <- length(x)
20          stopifnot(lx == 1 || lx == n) # but keep 'x' short for now          lx.1 <- lx == 1L
21            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
22          if(is.logical(x))          if(is.logical(x))
23              cl <- "ldiMatrix"              cl <- "ldiMatrix"
24          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 27  Line 28
28          else if(is.complex(x)) {          else if(is.complex(x)) {
29              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
30          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
31            if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
32                new(cl, Dim = c(n,n), diag = "U")
33            else
34          new(cl, Dim = c(n,n), diag = "N",          new(cl, Dim = c(n,n), diag = "N",
35              x = if(lx == 1) rep.int(x,n) else x)                  x = if(lx.1) rep.int(x,n) else x)
36      }      }
37  }  }
38
# Line 277  Line 281
281  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
282  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
283      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
284      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
285      y      y
286  }  }
287
# Line 294  Line 298
298                mod.x <- mode(x@x)                mod.x <- mode(x@x)
299                r <- vector(mod.x, length = n^2)                r <- vector(mod.x, length = n^2)
300                if(n)                if(n)
301                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1L) * (n + 1)] <-
302                        if(x@diag == "U") as1(mod=mod.x) else x@x                        if(x@diag == "U") as1(mod=mod.x) else x@x
303                r                r
304            })            })
# Line 869  Line 873
873  ##       Compare --> logical  ##       Compare --> logical
874  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
875
876  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
877  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
878  for(arg2 in c("numeric","logical"))  for(arg2 in c("numeric","logical"))
879  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
880            function(e1,e2) {            function(e1,e2) {
881                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
882                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
883                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
884                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
885                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
886                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
887                        e1@diag <- "N"                        e1@diag <- "N"
888                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
889                    } else                        } ## else: result = e1  (is "U" diag)
890                      } else {
891                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
892                    if(length(r))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
893                        e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
return(e1)
894                }                }
895                      e1
896                  } else
897                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
898            })            })
899
900  for(arg1 in c("numeric","logical"))  for(arg1 in c("numeric","logical"))
901  setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
902            function(e1,e2) {            function(e1,e2) {
903                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
904                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
905                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
906                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
907                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
908                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
909                        e2@diag <- "N"                        e2@diag <- "N"
910                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
911                    } else                        } ## else: result = e2  (is "U" diag)
912                      } else {
913                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
914                    if(length(r))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
915                        e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
return(e2)
916                }                }
917                      e2
918                  } else
919                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
920            })            })
921
# Line 915  Line 923
923  for(arg2 in c("numeric","logical"))  for(arg2 in c("numeric","logical"))
924  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
925            function(e1,e2) {            function(e1,e2) {
926                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
927                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
928                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
929                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
930                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
931                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
932                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
933                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
934                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
935                              E@diag <- "N"
936                              E@x[seq_len(n)] <- r # possibly recycling r
937                          } ## else: result = E  (is "U" diag)
938                      } else {
939                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
940                    ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix                    ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
941                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
if(length(r)) {
if(!is.double(r)) r <- as.double(r)
e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
}
return(e1)
942                }                }
943                      E
944                  } else
945                callGeneric(diag2tT.u(e1,e2, "l"), e2)                callGeneric(diag2tT.u(e1,e2, "l"), e2)
946            })            })
947
948  for(arg1 in c("numeric","logical"))  for(arg1 in c("numeric","logical"))
949  setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
950            function(e1,e2) {            function(e1,e2) {
951                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
952                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
953                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
954                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
955                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
956                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
957                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
958                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
959                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
960                              E@diag <- "N"
961                              E@x[seq_len(n)] <- r # possibly recycling r
962                          } ## else: result = E  (is "U" diag)
963                      } else {
964                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
965                    ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix                    ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
966                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
if(length(r)) {
if(!is.double(r)) r <- as.double(r)
e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
}
return(e2)
967                }                }
968                      E
969                  } else
970                callGeneric(e1, diag2tT.u(e2,e1, "l"))                callGeneric(e1, diag2tT.u(e2,e1, "l"))
971            })            })
972
973  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
974    ##
975    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
976    if(FALSE) {
977        selectMethod("<", c("numeric","lMatrix"))# Compare
978        selectMethod("&", c("numeric","lMatrix"))# Logic
979    } ## so we don't need to define a method here :
980
981  for(arg2 in c("numeric","logical"))  for(arg2 in c("numeric","logical"))
982  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
983            function(e1,e2) {            function(e1,e2) {
984                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
985                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
986                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
987                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
988                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
989                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
990                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
991                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
992                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
993                              E@diag <- "N"
994                              E@x[seq_len(n)] <- r # possibly recycling r
995                          } ## else: result = E  (is "U" diag)
996                      } else {
997                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
998                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
999                    if(length(r)) {                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
if(!is.logical(r)) r <- as.logical(r)
e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
}
return(e1)
1000                }                }
1001                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
})

for(arg1 in c("numeric","logical"))
setMethod("Ops", signature(e1 = arg1, e2 = "ddiMatrix"),
function(e1,e2) {
n <- e2@Dim[1]; nsq <- n^2
f0 <- callGeneric(e1, 0)
if(all(is0(f0))) { # remain diagonal
L1 <- (le <- length(e1)) == 1L
if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
e2@diag <- "N"
if(L1) r <- rep.int(r, n)
1002                    } else                    } else
1003                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "d"), e2)
if(length(r)) {
if(!is.logical(r)) r <- as.logical(r)
e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
}
e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
return(e2)
}
callGeneric(e1, diag2tT.u(e2,e1, "d"))
1004            })            })
1005
1006  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1007  for(arg2 in c("numeric","logical"))  for(arg2 in c("numeric","logical"))
1008  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1009            function(e1,e2) {            function(e1,e2) {
1010                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1011                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1012                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1013                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1014                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
1015                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1016                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1017                        e1@diag <- "N"                        e1@diag <- "N"
1018                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1019                    } else                        } ## else: result = e1  (is "U" diag)
1020                      } else {
1021                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1022                    if(length(r))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1023                        e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
return(e1)
1024                }                }
1025                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
})

for(arg1 in c("numeric","logical"))
setMethod("Ops", signature(e1 = arg1, e2 = "ldiMatrix"),
function(e1,e2) {
n <- e2@Dim[1]; nsq <- n^2
f0 <- callGeneric(e1, FALSE)
if(all(is0(f0))) { # remain diagonal
L1 <- (le <- length(e1)) == 1L
if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
e2@diag <- "N"
if(L1) r <- rep.int(r, n)
1026                    } else                    } else
1027                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "l"), e2)
if(length(r))
e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
return(e2)
}
callGeneric(e1, diag2tT.u(e2,e1, "l"))
1028            })            })
1029
1030

1031  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1032  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1033      ## ddi*:      ## ddi*:
# Line 1166  Line 1144
1144            })            })
1145
1146  rm(dense.subCl, diCls)# not used elsewhere  rm(dense.subCl, diCls)# not used elsewhere
1147
1148    setMethod("summary", signature(object = "diagonalMatrix"),
1149              function(object, ...) {
1150                  d <- dim(object)
1151                  r <- summary(object@x, ...)
1153                      sprintf('%d x %d diagonal Matrix of class "%s"',
1154                              d[1], d[2], class(object))
1155                  ## use ole' S3 technology for such a simple case
1156                  class(r) <- c("diagSummary", class(r))
1157                  r
1158              })
1159
1160    print.diagSummary <- function (x, ...) {