Diff of /partyMod/R/Utils.R [000000] .. [fbf06f]

Switch to side-by-side view

--- a
+++ b/partyMod/R/Utils.R
@@ -0,0 +1,248 @@
+
+# $Id$
+
+### Wrapper for functions defined in ./src/Utilsc
+
+qsvd <- function(x) {
+    if (!is.matrix(x) || ncol(x) != nrow(x))
+        stop(sQuote("x"), " is not a quadratic matrix")
+
+    svdmem <- new("svd_mem", ncol(x)) 
+    dummy <- .Call("R_svd", x, svdmem, PACKAGE = "atlantisPartyMod")
+    rm(dummy)
+    return(list(u = svdmem@u, vt = svdmem@v, d = svdmem@s))
+}
+
+MPinv <- function(x, tol = sqrt(.Machine$double.eps)) {
+    if (!is.matrix(x) || ncol(x) != nrow(x))
+        stop(sQuote("x"), " is not a quadratic matrix")
+
+    svdmem <- new("svd_mem", ncol(x))
+    RET <- .Call("R_MPinv", x, tol, svdmem, PACKAGE = "atlantisPartyMod")
+    return(RET@MPinv)
+}
+
+### wrapper for functions defined in ./src/Splits.c
+
+Split <- function(x, y, weights, splitctrl) {
+
+    if (is.factor(y))
+        ym <- sapply(levels(y), function(l) as.numeric(y == l))
+     else 
+        ym <- matrix(y, ncol = 1)
+    storage.mode(ym) <- "double"
+
+    if (is.factor(x)) {
+        xm <- sapply(levels(x), function(l) as.numeric(x == l))
+        storage.mode(xm) <- "double"
+        xc <- as.numeric(x)
+        storage.mode(xc) <- "integer"
+        lecxy <- new("LinStatExpectCovar", ncol(xm), ncol(ym))
+        lec <- new("LinStatExpectCovar", as.integer(1), ncol(ym))
+        eci <- ExpectCovarInfluence(ym, weights)
+        split <- .Call("R_splitcategorical", xm, xc, ym, weights, lec, lecxy,
+                       eci, splitctrl, PACKAGE = "atlantisPartyMod")
+    } else {
+        ox <- order(x)
+        storage.mode(ox) <- "integer"
+        xm <- matrix(x, ncol = 1)
+        storage.mode(xm) <- "double"
+        lec <- new("LinStatExpectCovar", as.integer(1), ncol(ym))
+        eci <- ExpectCovarInfluence(ym, weights)
+        split <- .Call("R_split", xm, ym, weights, ox, lec,
+                       eci, splitctrl, PACKAGE = "atlantisPartyMod")
+    }
+    split
+}
+
+### Wrapper for functions defined in ./src/TestStatistic.c
+
+maxabsTestStatistic <- function(t, mu, Sigma, tol = sqrt(.Machine$double.eps)) {
+    storage.mode(t) <- "double"
+    storage.mode(mu) <- "double"
+    storage.mode(Sigma) <- "double"
+    storage.mode(tol) <- "double"
+    
+    if (length(t) != length(mu) || length(t) != nrow(Sigma)) 
+        stop("dimensions don't match")
+        
+    .Call("R_maxabsTestStatistic", t, mu, Sigma, tol, PACKAGE = "atlantisPartyMod")
+}
+
+quadformTestStatistic <- function(t, mu, Sigma, 
+    tol = sqrt(.Machine$double.eps)) {
+
+    storage.mode(t) <- "double"
+    storage.mode(mu) <- "double"
+    storage.mode(Sigma) <- "double"
+    storage.mode(tol) <- "double"
+    
+    if (length(t) != length(mu) || length(t) != nrow(Sigma)) 
+        stop("dimensions don't match")
+        
+    SigmaPlus <- MPinv(Sigma, tol = tol)
+    .Call("R_quadformTestStatistic", t, mu, SigmaPlus, PACKAGE = "atlantisPartyMod")
+}
+
+
+### Wrapper for functions defined in ./src/LinearStatistic.c
+
+LinearStatistic <- function(x, y, weights) {
+    storage.mode(x) <- "double"
+    storage.mode(y) <- "double"
+    storage.mode(weights) <- "double"
+    .Call("R_LinearStatistic", x, y, weights, PACKAGE = "atlantisPartyMod")
+}
+
+ExpectCovarInfluence <- function(y, weights) {
+    storage.mode(y) <- "double"
+    storage.mode(weights) <- "double"
+    .Call("R_ExpectCovarInfluence", y, weights, PACKAGE = "atlantisPartyMod")
+}
+
+ExpectCovarLinearStatistic <- function(x, y, weights) {
+    storage.mode(x) <- "double"
+    storage.mode(y) <- "double"
+    storage.mode(weights) <- "double"
+    expcovinf <- ExpectCovarInfluence(y, weights)
+    .Call("R_ExpectCovarLinearStatistic", x, y, weights, expcovinf,
+          PACKAGE = "atlantisPartyMod")
+}
+
+PermutedLinearStatistic <- function(x, y, indx, perm) {
+    storage.mode(x) <- "double"
+    storage.mode(y) <- "double"
+
+    if (any(indx < 1 || indx > nrow(y)))
+        stop("wrong indices")
+
+    if (any(perm < 1 || perm > nrow(y)))
+        stop("wrong indices")
+
+    # C indexing 
+    indx <- indx - 1
+    perm <- perm - 1
+    storage.mode(indx) <- "integer"
+    storage.mode(perm) <- "integer"
+    .Call("R_PermutedLinearStatistic", x, y, indx, perm, 
+          PACKAGE = "atlantisPartyMod")
+}
+
+### Median Survival Time, see survival:::print.survfit
+
+mst <- function(x) {
+    minmin <- function(y, xx) {
+
+        ### don't complain about Inf 
+        ww <- getOption("warn")
+        on.exit(options(warn = ww))
+        options(warn = -1)
+
+        if (any(!is.na(y) & y==.5)) {
+            if (any(!is.na(y) & y <.5))
+                .5*(min(xx[!is.na(y) & y==.5]) + min(xx[!is.na(y) & y<.5]))
+            else
+                .5*(min(xx[!is.na(y) & y==.5]) + max(xx[!is.na(y) & y==.5]))
+        } else   min(xx[!is.na(y) & y<=.5])
+    }
+    med <- minmin(x$surv, x$time)
+    return(med)
+}
+
+dostep <- function(x, y) {
+
+    ### create a step function based on x, y coordinates
+    ### modified from `survival:print.survfit'
+    if (is.na(x[1] + y[1])) {
+        x <- x[-1]
+        y <- y[-1]
+    }
+    n <- length(x)
+    if (n > 2) {  
+        # replace verbose horizonal sequences like
+        # (1, .2), (1.4, .2), (1.8, .2), (2.3, .2), (2.9, .2), (3, .1)
+        # with (1, .2), (3, .1).  They are slow, and can smear the looks
+        # of the line type.
+        dupy <- c(TRUE, diff(y[-n]) !=0, TRUE)
+        n2 <- sum(dupy)
+
+        #create a step function
+        xrep <- rep(x[dupy], c(1, rep(2, n2-1)))
+        yrep <- rep(y[dupy], c(rep(2, n2-1), 1))
+        RET <- list(x = xrep, y = yrep)
+    } else {
+        if (n == 1) {
+            RET <- list(x = x, y = y)
+        } else {
+            RET <- list(x = x[c(1,2,2)], y = y[c(1,1,2)])
+        }
+    }
+    return(RET)
+}
+
+### Normalized mutual information, Stehl & Gosh (2002) JMLR
+nmi <- function(x, y)
+{
+    x <- table(x, y)
+    x <- x / sum(x)                     # ???
+
+    m_x <- rowSums(x)
+    m_y <- colSums(x)
+    y <- outer(m_x, m_y)
+
+    i <- which((x > 0) & (y > 0))
+    out <- sum(x[i] * log(x[i] / y[i]))
+    e_x <- sum(m_x * log(ifelse(m_x > 0, m_x, 1)))
+    e_y <- sum(m_y * log(ifelse(m_y > 0, m_y, 1)))
+
+    out / sqrt(e_x * e_y)
+}
+
+### check if two objects are identical and print differences else
+isequal <- function(a, b) {
+    attributes(a) <- NULL
+    attributes(b) <- NULL
+    if (!isTRUE(all.equal(a, b))) {
+        print(a, digits = 10)
+        print(b, digits = 10)
+        return(FALSE)
+    } else {
+        return(TRUE)
+    }
+}
+
+mysurvfit <- function(y, weights, ...) {
+
+    stopifnot(extends(class(y), "Surv"))
+    ### see comment on weights and subset in ?survfit
+    y <- y[weights > 0,]
+    weights <- weights[weights > 0]
+    return(survfit(y ~ 1, weights = weights, ...))
+}
+
+R_get_nodeID <- function(tree, inputs, mincriterion)
+    .Call("R_get_nodeID", tree, inputs, 0.0, -1L, PACKAGE = "atlantisPartyMod")
+
+R_getpredictions <- function(tree, where)
+    .Call("R_getpredictions", tree, where, PACKAGE = "atlantisPartyMod")
+
+R_remove_weights <- function(tree, remove)
+    .Call("R_remove_weights", tree, remove, package = "atlantisPartyMod")
+
+R_modify_response <- function(y, responses)
+    .Call("R_modify_response", as.double(y), responses,
+          PACKAGE = "atlantisPartyMod")
+
+R_TreeGrow <- function(object, weights, fitmem, ctrl, where)
+    .Call("R_TreeGrow", object, weights, fitmem, ctrl,
+          where, PACKAGE = "atlantisPartyMod")
+
+copyslots <- function(source, target) {
+    slots <- names(getSlots(class(source)))
+    slots <- slots[(slots %in% names(getSlots(class(target))))]
+    if (length(slots) == 0) 
+        stop("no common slots to copy to")
+    for (s in slots)
+        eval(parse(text = paste("target@", s, " <- source@", s)))
+    return(target)
+}