SCM

SCM Repository

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

Diff of /pkg/R/Auxiliaries.R

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

revision 703, Thu Apr 21 22:00:04 2005 UTC revision 1654, Fri Oct 27 16:58:15 2006 UTC
# Line 1  Line 1 
1  #### "Namespace private" Auxiliaries  such as method functions  #### "Namespace private" Auxiliaries  such as method functions
2  #### (called from more than one place --> need to be defined early)  #### (called from more than one place --> need to be defined early)
3    
4    .isR_24 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.4")
5    
6    ## Need to consider NAs ;  "== 0" even works for logical & complex:
7    is0  <- function(x) !is.na(x) & x == 0
8    isN0 <- function(x)  is.na(x) | x != 0
9    all0 <- function(x) !any(is.na(x)) && all(x == 0)
10    
11    allTrue  <- function(x) !any(is.na(x)) && all(x)
12    allFalse <- function(x) !any(is.na(x)) && !any(x)
13    
14    
15  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):
16  .M.v <- function(x, y) callGeneric(x, as.matrix(y))  .M.v <- function(x, y) callGeneric(x, as.matrix(y))
17  .v.M <- function(x, y) callGeneric(rbind(x), y)  .v.M <- function(x, y) callGeneric(rbind(x), y)
18    
19    .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)
20    
21  .has.DN <- ## has non-trivial Dimnames slot?  .has.DN <- ## has non-trivial Dimnames slot?
22      function(x) !identical(list(NULL,NULL), x@Dimnames)      function(x) !identical(list(NULL,NULL), x@Dimnames)
23    
24    .bail.out.1 <- function(fun, cl) {
25        stop(gettextf('not-yet-implemented method for %s(<%s>).\n ->>  Ask the package authors to implement the missing feature.', fun, cl),
26             call. = FALSE)
27    }
28    .bail.out.2 <- function(fun, cl1, cl2) {
29        stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>).\n ->>  Ask the package authors to implement the missing feature.',
30                      fun, cl1, cl2), call. = FALSE)
31    }
32    
33    ## This should be done in C and be exported by 'methods':  [FIXME - ask JMC ]
34    copyClass <- function(x, newCl, sNames =
35                          intersect(slotNames(newCl), slotNames(x))) {
36        r <- new(newCl)
37        for(n in sNames)
38            slot(r, n) <- slot(x, n)
39        r
40    }
41    
42  ## chol() via "dpoMatrix"  ## chol() via "dpoMatrix"
43  cholMat <- function(x, pivot, LINPACK) {  cholMat <- function(x, pivot, ...) {
44      px <- as(x, "dpoMatrix")      px <- as(x, "dpoMatrix")
45      if (isTRUE(validObject(px, test=TRUE))) chol(px)      if (isTRUE(validObject(px, test=TRUE))) chol(px)
46      else stop("'x' is not positive definite -- chol() undefined.")      else stop("'x' is not positive definite -- chol() undefined.")
47  }  }
48    
49    dimCheck <- function(a, b) {
50        da <- dim(a)
51        db <- dim(b)
52        if(any(da != db))
53            stop(gettextf("Matrices must have same dimensions in %s",
54                          deparse(sys.call(sys.parent()))),
55                 call. = FALSE)
56        da
57    }
58    
59    dimNamesCheck <- function(a, b) {
60        ## assume dimCheck() has happened before
61        nullDN <- list(NULL,NULL)
62        h.a <- !identical(nullDN, dna <- dimnames(a))
63        h.b <- !identical(nullDN, dnb <- dimnames(b))
64        if(h.a || h.b) {
65            if (!h.b) dna
66            else if(!h.a) dnb
67            else { ## both have non-trivial dimnames
68                r <- dna # "default" result
69                for(j in 1:2) {
70                    dn <- dnb[[j]]
71                    if(is.null(r[[j]]))
72                        r[[j]] <- dn
73                    else if (!is.null(dn) && any(r[[j]] != dn))
74                        warning(gettextf("dimnames [%d] mismatch in %s", j,
75                                         deparse(sys.call(sys.parent()))),
76                                call. = FALSE)
77                }
78                r
79            }
80        }
81        else
82            nullDN
83    }
84    
85    rowCheck <- function(a, b) {
86        da <- dim(a)
87        db <- dim(b)
88        if(da[1] != db[1])
89            stop(gettextf("Matrices must have same number of rows in %s",
90                          deparse(sys.call(sys.parent()))),
91                 call. = FALSE)
92        ## return the common nrow()
93        da[1]
94    }
95    
96    colCheck <- function(a, b) {
97        da <- dim(a)
98        db <- dim(b)
99        if(da[2] != db[2])
100            stop(gettextf("Matrices must have same number of columns in %s",
101                          deparse(sys.call(sys.parent()))),
102                 call. = FALSE)
103        ## return the common ncol()
104        da[2]
105    }
106    
107    ## Note: !isPacked(.)  i.e. `full' still contains
108    ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
109    isPacked <- function(x)
110    {
111        ## Is 'x' a packed (dense) matrix ?
112        is(x, "denseMatrix") &&
113        any("x" == slotNames(x)) && length(x@x) < prod(dim(x))
114    }
115    
116    emptyColnames <- function(x)
117    {
118        ## Useful for compact printing of (parts) of sparse matrices
119        ## possibly  dimnames(x) "==" NULL :
120        dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
121        x
122    }
123    
124    ### TODO:  write in C and port to base (or 'utils') R
125    indTri <- function(n, upper = TRUE) {
126        ## == which(upper.tri(diag(n)) or
127        ##    which(lower.tri(diag(n)) -- but much more efficiently for largish 'n'
128        stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0)
129        if(n <= 2)
130            return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0))
131        ## First, compute the 'diff(.)'  fast.  Use integers
132        one <- 1:1 ; two <- 2:2
133        n1 <- n - one
134        n2 <- n1 - one
135        r <- rep.int(one, n*n1/two - one)
136        r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n
137        ## now have "dliu" difference; revert to "liu":
138        cumsum(c(if(upper) n+one else two, r))
139    }
140    
141    
142    prTriang <- function(x, digits = getOption("digits"),
143                         maxp = getOption("max.print"),
144                         justify = "none", right = TRUE)
145    {
146        ## modeled along stats:::print.dist
147        upper <- x@uplo == "U"
148    
149        m <- as(x, "matrix")
150        cf <- format(m, digits = digits, justify = justify)
151        if(upper)
152            cf[row(cf) > col(cf)] <- "."
153        else
154            cf[row(cf) < col(cf)] <- "."
155        if(.isR_24)
156             print(cf, quote = FALSE, right = right, max = maxp)
157        else print(cf, quote = FALSE, right = right)
158        invisible(x)
159    }
160    
161    prMatrix <- function(x, digits = getOption("digits"),
162                         maxp = getOption("max.print")) {
163        d <- dim(x)
164        cl <- class(x)
165        cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
166        if(prod(d) <= maxp) {
167            if(is(x, "triangularMatrix"))
168                prTriang(x, digits = digits, maxp = maxp)
169            else {
170                if(.isR_24)
171                     print(as(x, "matrix"), digits = digits, max = maxp)
172                else print(as(x, "matrix"), digits = digits)
173            }
174        }
175        else { ## d[1] > maxp / d[2] >= nr :
176            m <- as(x, "matrix")
177            nr <- maxp %/% d[2]
178            n2 <- ceiling(nr / 2)
179            print(head(m, max(1, n2)))
180            cat("\n ..........\n\n")
181            print(tail(m, max(1, nr - n2)))
182        }
183        ## DEBUG: cat("str(.):\n") ; str(x)
184        invisible(x)# as print() S3 methods do
185    }
186    
187    nonFALSE <- function(x) {
188        ## typically used for lMatrices:  (TRUE,NA,FALSE) |-> (TRUE,FALSE)
189        if(any(ix <- is.na(x))) x[ix] <- TRUE
190        x
191    }
192    
193    nz.NA <- function(x, na.value) {
194        ## Non-Zeros of x
195        ## na.value: TRUE: NA's give TRUE, they are not 0
196        ##             NA: NA's are not known ==> result := NA
197        ##          FALSE: NA's give FALSE, could be 0
198        stopifnot(is.logical(na.value) && length(na.value) == 1)
199        if(is.na(na.value)) x != 0
200        else  if(na.value)  isN0(x)
201        else                x != 0 & !is.na(x)
202    }
203    
204    ## Number of non-zeros :
205    ## FIXME? -- make this into a generic function (?)
206    nnzero <- function(x, na.counted = NA) {
207        ## na.counted: TRUE: NA's are counted, they are not 0
208        ##               NA: NA's are not known (0 or not) ==>  result := NA
209        ##            FALSE: NA's are omitted before counting
210        cl <- class(x)
211        if(!extends(cl, "Matrix"))
212            sum(nz.NA(x, na.counted))
213        else if(extends(cl, "sparseMatrix"))
214            ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
215           switch(.sp.class(cl),
216                   "CsparseMatrix" = length(x@i),
217                   "TsparseMatrix" = length(x@i),
218                   "RsparseMatrix" = length(x@j))
219        else ## denseMatrix
220            sum(nz.NA(as_geClass(x, cl)@x, na.counted))
221    }
222    
223    ## For sparseness handling
224    ## return a 2-column (i,j) matrix of
225    ## 0-based indices of non-zero entries  :
226    non0ind <- function(x) {
227    
228        if(is.numeric(x))
229            return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
230        ## else
231        stopifnot(is(x, "sparseMatrix"))
232        non0.i <- function(M) {
233            if(is(M, "TsparseMatrix"))
234                return(unique(cbind(M@i,M@j)))
235            if(is(M, "pMatrix"))
236                return(cbind(seq_len(nrow(M)), M@perm) - 1:1)
237            ## else:
238            isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
239            .Call(compressed_non_0_ij, M, isC)
240        }
241    
242        if(is(x, "symmetricMatrix")) { # also get "other" triangle
243            ij <- non0.i(x)
244            notdiag <- ij[,1] != ij[,2]# but not the diagonals again
245            rbind(ij, ij[notdiag, 2:1])
246        }
247        else if(is(x, "triangularMatrix")) { # check for "U" diag
248            if(x@diag == "U") {
249                i <- seq_len(dim(x)[1]) - 1:1
250                rbind(non0.i(x), cbind(i,i))
251            } else non0.i(x)
252        }
253        else
254            non0.i(x)
255    }
256    
257    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
258    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
259    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
260    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
261    
262    complementInd <- function(ij, dim)
263    {
264        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
265        ##          but as 1-based indices
266        n <- prod(dim)
267        if(n == 0) return(integer(0))
268        ii <- 1:n
269        ii[-(1 + encodeInd(ij, nr = dim[1]))]
270    }
271    
272    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
273    
274    intersectInd <- function(ij1, ij2, nrow) {
275        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
276        ## return only the *common* entries
277        decodeInd(intersect(encodeInd(ij1, nrow),
278                            encodeInd(ij2, nrow)), nrow)
279    }
280    
281    WhichintersectInd <- function(ij1, ij2, nrow) {
282        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
283        ## find *where*  common entries are in ij1 & ij2
284        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
285        ni <- !is.na(m1)
286        list(which(ni), m1[ni])
287    }
288    
289    
290    ### There is a test on this in ../tests/dgTMatrix.R !
291    
292    uniqTsparse <- function(x, class.x = c(class(x))) {
293        ## Purpose: produce a *unique* triplet representation:
294        ##          by having (i,j) sorted and unique
295        ## -----------------------------------------------------------
296        ## The following is not quite efficient {but easy to program,
297        ## and as() are based on C code  (all of them?)
298        ##
299        ## FIXME: Do it fast for the case where 'x' is already 'uniq'
300    
301        switch(class.x,
302               "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
303               "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
304               "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
305               ## do we need this for "logical" ones, there's no sum() there!
306               "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
307               "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
308               "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
309               ## do we need this for "logical" ones, there's no sum() there!
310               "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"),
311               "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"),
312               "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"),
313               ## otherwise:
314               stop("not yet implemented for class ", class.x))
315    }
316    
317    ## Note: maybe, using
318    ## ----    xj <- .Call(Matrix_expand_pointers, x@p)
319    ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
320    ## but really efficient would be to use only one .Call(.) for uniq(.) !
321    
322    drop0 <- function(x, clx = c(class(x))) {
323        ## FIXME: Csparse_drop should do this (not losing symm./triang.):
324        ## Careful: 'Csparse_drop' also drops triangularity,...
325        ## .Call(Csparse_drop, as_CspClass(x, clx), 0)
326        as_CspClass(.Call(Csparse_drop, as_CspClass(x, clx), 0.),
327                    clx)
328    }
329    
330    uniq <- function(x) {
331        if(is(x, "TsparseMatrix")) uniqTsparse(x) else
332        if(is(x, "sparseMatrix")) drop0(x) else x
333    }
334    
335    asTuniq <- function(x) {
336        if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
337    }
338    
339    if(FALSE) ## try an "efficient" version
340    uniq_gT <- function(x)
341    {
342        ## Purpose: produce a *unique* triplet representation:
343        ##          by having (i,j) sorted and unique
344        ## ------------------------------------------------------------------
345        ## Arguments: a "gT" Matrix
346        stopifnot(is(x, "gTMatrix"))
347        if((n <- length(x@i)) == 0) return(x)
348        ii <- order(x@i, x@j)
349        if(any(ii != 1:n)) {
350            x@i <- x@i[ii]
351            x@j <- x@j[ii]
352            x@x <- x@x[ii]
353        }
354        ij <- x@i + nrow(x) * x@j
355        if(any(dup <- duplicated(ij))) {
356    
357        }
358        ### We should use a .Call() based utility for this!
359    
360    }
361    
362    t_geMatrix <- function(x) {
363        x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
364        x@Dim <- x@Dim[2:1]
365        x@Dimnames <- x@Dimnames[2:1]
366        ## FIXME: how to set factors?
367        x
368    }
369    
370    ## t( [dl]trMatrix ) and  t( [dl]syMatrix ) :
371    t_trMatrix <- function(x) {
372        x@x <- as.vector(t(as(x, "matrix")))
373        x@Dim <- x@Dim[2:1]
374        x@Dimnames <- x@Dimnames[2:1]
375        x@uplo <- if (x@uplo == "U") "L" else "U"
376        # and keep x@diag
377        x
378    }
379    
380    fixupDense <- function(m, from) {
381        if(is(m, "triangularMatrix")) {
382            m@uplo <- from@uplo
383            m@diag <- from@diag
384        } else if(is(m, "symmetricMatrix")) {
385            m@uplo <- from@uplo
386        }
387        m
388    }
389    
390    ## -> ./ldenseMatrix.R :
391    l2d_Matrix <- function(from) {
392        stopifnot(is(from, "lMatrix"))
393        fixupDense(new(sub("^l", "d", class(from)),
394                       x = as.double(from@x),
395                       Dim = from@Dim, Dimnames = from@Dimnames),
396                   from)
397        ## FIXME: treat 'factors' smartly {not for triangular!}
398    }
399    
400    ## -> ./ndenseMatrix.R :
401    n2d_Matrix <- function(from) {
402        stopifnot(is(from, "nMatrix"))
403        fixupDense(new(sub("^n", "d", class(from)),
404                       x = as.double(from@x),
405                       Dim = from@Dim, Dimnames = from@Dimnames),
406                   from)
407        ## FIXME: treat 'factors' smartly {not for triangular!}
408    }
409    n2l_spMatrix <- function(from) {
410        stopifnot(is(from, "nMatrix"))
411        new(sub("^n", "l", class(from)),
412            ##x = as.double(from@x),
413            Dim = from@Dim, Dimnames = from@Dimnames)
414    }
415    
416    if(FALSE)# unused
417    l2d_meth <- function(x) {
418        cl <- class(x)
419        as(callGeneric(as(x, sub("^l", "d", cl))), cl)
420    }
421    
422    ## return "d" or "l" or "n" or "z"
423    .M.kind <- function(x, clx = class(x)) {
424        if(is.matrix(x)) { ## 'old style matrix'
425            if     (is.numeric(x)) "d"
426            else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
427            else if(is.complex(x)) "z"
428            else stop("not yet implemented for matrix w/ typeof ", typeof(x))
429        }
430        else if(extends(clx, "dMatrix")) "d"
431        else if(extends(clx, "nMatrix")) "n"
432        else if(extends(clx, "lMatrix")) "l"
433        else if(extends(clx, "zMatrix")) "z"
434        else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
435        else stop(" not yet be implemented for ", clx)
436    }
437    
438    .M.shape <- function(x, clx = class(x)) {
439        if(is.matrix(x)) { ## 'old style matrix'
440            if     (isDiagonal  (x)) "d"
441            else if(isTriangular(x)) "t"
442            else if(isSymmetric (x)) "s"
443            else "g" # general
444        }
445        else if(extends(clx, "diagonalMatrix"))  "d"
446        else if(extends(clx, "triangularMatrix"))"t"
447        else if(extends(clx, "symmetricMatrix")) "s"
448        else "g"
449    }
450    
451    
452    class2 <- function(cl, kind = "l", do.sub = TRUE) {
453        ## Find "corresponding" class; since pos.def. matrices have no pendant:
454        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
455        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
456        else if(do.sub) sub("^d", kind, cl)
457        else cl
458    }
459    
460    geClass <- function(x) {
461        if     (is(x, "dMatrix")) "dgeMatrix"
462        else if(is(x, "lMatrix")) "lgeMatrix"
463        else if(is(x, "nMatrix")) "ngeMatrix"
464        else if(is(x, "zMatrix")) "zgeMatrix"
465        else stop("general Matrix class not yet implemented for ",
466                  class(x))
467    }
468    
469    .dense.prefixes <- c("d" = "di",
470                         "t" = "tr",
471                         "s" = "sy",
472                         "g" = "ge")
473    
474    .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
475                          "t" = "t",
476                          "s" = "s",
477                          "g" = "g")
478    
479    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
480    as_dense <- function(x) {
481        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
482    }
483    
484    .sp.class <- function(x) { ## find and return the "sparseness class"
485        if(!is.character(x)) x <- class(x)
486        for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
487            if(extends(x, cl))
488                return(cl)
489        ## else (should rarely happen)
490        as.character(NA)
491    }
492    
493    as_Csparse <- function(x) {
494        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
495    }
496    as_Rsparse <- function(x) {
497        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
498    }
499    as_Tsparse <- function(x) {
500        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
501    }
502    
503    as_geClass <- function(x, cl) {
504        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
505            as(x, cl)
506        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x)) {
507            kind <- .M.kind(x)
508            as(x, class2(cl, kind, do.sub= kind != "d"))
509        } else if(extends(cl, "triangularMatrix") && isTriangular(x))
510            as(x, cl)
511        else
512            as(x, paste(.M.kind(x), "geMatrix", sep=''))
513    }
514    
515    as_CspClass <- function(x, cl) {
516        if (## diagonal is *not* sparse:
517            ##(extends(cl, "diagonalMatrix") && isDiagonal(x)) ||
518            (extends(cl, "symmetricMatrix") && isSymmetric(x)) ||
519            (extends(cl, "triangularMatrix")&& isTriangular(x)))
520            as(x, cl)
521        else if(is(x, "CsparseMatrix")) x
522        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
523    }
524    
525    
526    ## -> ./ddenseMatrix.R :
527    d2l_Matrix <- function(from) {
528        stopifnot(is(from, "dMatrix"))
529        fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
530                       Dim = from@Dim, Dimnames = from@Dimnames),
531                   from)
532        ## FIXME: treat 'factors' smartly {not for triangular!}
533    }
534    
535    
536    try_as <- function(x, classes, tryAnyway = FALSE) {
537        if(!tryAnyway && !is(x, "Matrix"))
538            return(x)
539        ## else
540        ok <- canCoerce(x, classes[1])
541        while(!ok && length(classes <- classes[-1])) {
542            ok <- canCoerce(x, classes[1])
543        }
544        if(ok) as(x, classes[1]) else x
545    }
546    
547    
548    ## For *dense* matrices
549    isTriMat <- function(object, upper = NA) {
550        ## pretest: is it square?
551        d <- dim(object)
552        if(d[1] != d[2]) return(FALSE)
553        ## else slower test
554        if(!is.matrix(object))
555            object <- as(object,"matrix")
556        if(is.na(upper)) {
557            if(all0(object[lower.tri(object)]))
558                structure(TRUE, kind = "U")
559            else if(all0(object[upper.tri(object)]))
560                structure(TRUE, kind = "L")
561            else FALSE
562        } else if(upper)
563            all0(object[lower.tri(object)])
564        else ## upper is FALSE
565            all0(object[upper.tri(object)])
566    }
567    
568    ## For Csparse matrices
569    isTriC <- function(x, upper = NA) {
570        ## pretest: is it square?
571        d <- dim(x)
572        if(d[1] != d[2]) return(FALSE)
573        ## else
574        if(d[1] == 0) return(TRUE)
575        ni <- 1:d[2]
576        ## the row indices split according to column:
577        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
578        lil <- unlist(lapply(ilist, length), use.names = FALSE)
579        if(any(lil == 0)) {
580            pos <- lil > 0
581            if(!any(pos)) ## matrix of all 0's
582                return(TRUE)
583            ilist <- ilist[pos]
584            ni <- ni[pos]
585        }
586        ni0 <- ni - 1:1 # '0-based ni'
587        if(is.na(upper)) {
588            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
589                structure(TRUE, kind = "U")
590            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
591                structure(TRUE, kind = "L")
592            else FALSE
593        } else if(upper) {
594            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
595        } else { ## 'lower'
596            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
597        }
598    }
599    
600    .is.diagonal <- function(object) {
601        ## "matrix" or "denseMatrix" (but not "diagonalMatrix")
602        d <- dim(object)
603        if(d[1] != (n <- d[2])) FALSE
604        else if(is.matrix(object))
605            ## requires that "vector-indexing" works for 'object' :
606            all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
607        else ## "denseMatrix" -- packed or unpacked
608            if(is(object, "generalMatrix")) # "dge", "lge", ...
609                all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
610            else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
611                ## -> has 'uplo'  differentiate between packed and unpacked
612    
613    ### .......... FIXME ...............
614    
615                packed <- isPacked(object)
616                if(object@uplo == "U") {
617                } else { ## uplo == "L"
618                }
619    
620    ### very cheap workaround
621                all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
622            }
623    }
624    
625    
626    ## FIXME? -- this should also work for "ltT", "ntT", ... :
627    diagU2N <- function(x)
628    {
629        ## Purpose: Transform a *unit diagonal* sparse triangular matrix
630        ##  into one with explicit diagonal entries '1'
631        xT <- as(x, "dgTMatrix")
632        ## leave it as  T* - the caller can always coerce to C* if needed:
633        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
634            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
635    }
636    
637    ## FIXME: this should probably be dropped / replaced by as_Csparse
638    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
639        x <- as(x, "dgCMatrix")
640        callGeneric()
641    }
642    
643    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
644        ## used e.g. inside colSums() etc methods
645        x <- as(x, "dgTMatrix")
646        callGeneric()
647    }
648    
649    
650    ### Fast much simplified version of tapply()
651    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
652        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
653    }
654    
655    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
656    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
657    ## }
658    

Legend:
Removed from v.703  
changed lines
  Added in v.1654

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