R/clust.util.R

In immunoClust: immunoClust - Automated Pipeline for Population Detection in Flow Cytometry

###
##
bhattacharyya.prob <- function(gM,gS, cM,cS, alpha=1)
{
    ret <- 0
    if( alpha <= 0 ) {
        ret <- bhattacharyya.prob(gM,diag(diag(gS)), cM, diag(diag(cS)))
        return(ret)
    }
    else if( alpha < 1 ) {
        a <- bhattacharyya.prob(gM,gS, cM, cS)
        b <- bhattacharyya.prob(gM,diag(diag(gS)), cM, diag(diag(cS)))
        ret <- alpha*a + (1-alpha)*b
        return(ret)
    }
    
    det_g <- log(det(gS))
    det_c <- log(det(cS))
    
    S <- chol((gS+cS)/2)
    S <- chol2inv(S)
    
    logD <- log(det(S))
    dist <- mahalanobis(gM, cM, S, inverted=TRUE)
    logD <- logD + 0.5*det_g + 0.5*det_c
    logD <- logD - 0.5*0.5*dist
    
    ## normalization factor
    logD <- logD - 0.25*det_g
    
    ret <- exp(0.5*logD)
    if( is.na(ret) || ret <= 0 ) {
        ret <- 0
    }
    
    ret
    
}

## bhattacharyya.dist = -log(bhattacharyya.coeff)
bhattacharyya.dist <- function (gM, gS, cM, cS)
{
    if( is.null(gS) || is.null(cS) ) {
        
        return(0)
    }
    
    S <- (gS + cS)/2
    d1 <- mahalanobis(gM, cM, S)/8
    
    d2 <- log(det(as.matrix(S))/sqrt(det(as.matrix(gS)) *
    det(as.matrix(cS))))/2
    ret <- d1 + d2
    ret
}

## bhattacharyya.coeff = exp(-bhattacharyya.dist)
bhattacharyya.coeff <- function(gM,gS, cM,cS, alpha=1)
{
    ret <- 0
    if( is.null(gS) || is.null(cS) ) {
        
        return(0)
    }
    if( alpha <= 0 ) {
        ret <- bhattacharyya.coeff(gM,diag(diag(gS)), cM, diag(diag(cS)))
        return(ret)
    }
    else if( alpha < 1 ) {
        a <- bhattacharyya.coeff(gM,gS, cM, cS)
        b <- bhattacharyya.coeff(gM,diag(diag(gS)), cM, diag(diag(cS)))
        ret <- alpha*a + (1-alpha)*b
        return(ret)
    }
    
    
    det_g <- log(det(gS))
    det_c <- log(det(cS))
    
    S <- NULL
    try( S <- solve((gS+cS) / 2), silent=TRUE)
    
    if( is.null(S) ) {
        return (0)
    }
    
    logD <- log(det(S))
    dist <- mahalanobis(gM, cM, S, inverted=TRUE)
    logD <- logD + 0.5*det_g + 0.5*det_c
    logD <- logD - 0.5*0.5*dist
    
    ret <- exp(0.5*logD)
    if( is.na(ret) || ret <= 0 ) {
        ret <- 0
    }
    
    if( ret > 1 ) {
        ret <- 1
    }
    
    ret
    
}
####

###
## clust.hclass
##  retrieve group membership from hcPairs structure
###
.clust.partconv <- function(x, consec = TRUE)
{
    n <- length(x)
    y <- numeric(n)
    u <- unique(x)
    if(consec) {
# number groups in order of first row appearance
        l <- length(u)
        for(i in seq_len(l))
        y[x == u[i]] <- i
    }
    else {
# represent each group by its lowest-numbered member
        for(i in u) {
            l <- x == i
            y[l] <- (1:n)[l][1]
        }
    }
    y
}
.clust.hclass <- function (hcPairs, G) 
{
    initial <- attributes(hcPairs)$init
    n <- length(initial)
    k <- length(unique(initial))
    G <- if (missing(G)) k:2 else rev(sort(unique(G)))
    select <- k - G
    
    if (length(select) == 1 && !select) 
    return(matrix(initial, ncol = 1, dimnames = list(NULL, as.character(G))))
    
    bad <- select < 0 | select >= k
    if (any(bad)) 
    stop("Some specified number of clusters are inconsistent")
    
## number of groupings to be labeled
    L <- length(select)
    cl <- matrix(NA, nrow = n, ncol = L, dimnames = list(NULL, as.character(G)))
    if (select[1]) { 
        m <- 1
    }
    else {
## highest G == n
        cl[, 1] <- initial
        m <- 2
    }
    for ( l in seq_len(max(select)) ) {
# merge at l    
        ij <- hcPairs[, l]
        i <- ij[1]
        j <- ij[2]
# i < j: all j became i
        initial[initial == j] <- i
        if (select[m] == l) {
            cl[, m] <- initial
            m <- m + 1
        }
    }
    apply(cl[, L:1, drop = FALSE], 2, .clust.partconv, consec = TRUE)
}
## clust.hclass

###
## convert hcPairs structure to hclust class
#.clust.hclust2 <- function(hcPairs)
#{
# li <- hcPairs[1,]
# lj <- hcPairs[2,]
# crit <- attributes(hcPairs)$change
# n <- length(li)
# obj <- .Fortran("hcass2", n=as.integer(n+1), 
#         ia=as.integer(li), ib=as.integer(lj),
#         order = integer(n+1), iia = integer(n+1), iib = integer(n+1),
#         PACKAGE="stats")
# hcls <- list(merge = cbind(obj$iia[1L:n], obj$iib[1L:n]), 
#        height = crit, crit=crit, order = obj$order, 
#        labels = 1:n)
# class(hcls) <- "hclust"
# hcls
#}
###

###
## convert hclust merging nodes to hcPairs merging nodes 
.clust.hpairs2 <- function(li, lj)
{
    hcp <- cbind(li,lj)
    N <- nrow(hcp)
    if( N>1)
    for( n in seq_len(N-1) ) {
        i <- hcp[n,1]
        j <- hcp[n,2]
        if( i > j ) {
            hcp[n,1] <- j
            hcp[n,2] <- i
            k <- j
        }
        else {
            k <- i
        }
        
        for( l in (n+1):N ) {
            if( hcp[l,1] == -n ) {
                hcp[l,1] <- k 
            }
            if( hcp[l,2] == -n ) {
                hcp[l,2] <- k
            }
        }
        
    }
    hcp
}
###

###
## 
.clust.mergedClusters <- function(x, cls)
{
    P <- ncol(x@mu)
    M <- rep(0, P)
    S <- matrix(0, nrow=P, ncol=P)
    
    if( length(cls) > 1 ) {
        M <- colSums(x@w[cls] * x@mu[cls,]) / sum(x@w[cls])
    }
    else {
        M <- x@mu[cls,]
    }
    
    for( i in cls ) {
        for( p in seq_len(P) ) {
            for( q in seq_len(P) ) {
                S[p,q] <- S[p,q] + (x@w[i]) * ( x@sigma[i, p, q] + 
                                    (x@mu[i,p]-M[p])*(x@mu[i,q]-M[q]) ) 
            }
        }
    }    
    S <- S/sum(x@w[cls])
    
    list("mu"=M,"sigma"=S)
}

###
##
.immunoClust2 <- function(obj, K, P, N, state=NULL,
expName="", parameters=c(), inc=c())
{
    L <- obj$L
    
#output sigma    
    sigma <- array(0, c(L, P, P))
    s <- matrix(obj$s, K, P * P, byrow=TRUE)
    for (k in seq_len(L))
    sigma[k,,] <- matrix(s[k,], P, P, byrow = TRUE)

    
# output mu    
    mu <- matrix(obj$m, K, P, byrow=TRUE)[seq_len(L),]
    dim(mu) <- c(L,P)
    
    
# output z
    if( length(inc) > 0 ) {
        z <- matrix(NA, length(inc), K)
        z[inc,] <- matrix(obj$z, N, K, byrow=TRUE)
        z <- z[,seq_len(L)]
        
        label <- rep(NA, length(inc))
        label[inc] <- obj$label
        
        N <- length(inc)
        
    }
    else {
        z <- matrix(obj$z, N, K, byrow=TRUE)
        z <- as.matrix(z[,seq_len(L)])
## 2014.07.08: quick fix (na/nan should not occur)  
        if( sum(is.infinite(z) | is.na(z) | is.nan(z)) > 0 ) {
            warning("Fehler: Z has infinite values", 
                sum(is.infinite(z) | is.na(z) | is.nan(z) ), "\n")
            z[is.infinite(z) | is.na(z) | is.nan(z)] <- 0
        }
        
        label <- obj$label
    }
    
    dim(z) <- c(N,L)
    
    
    if( !is.null(state) ) {
        for( k in seq_len(L) ) {
            state[k] <- state[ obj$history[k] ]
        }
        state <- state[seq_len(L)]
    }
    else {
        state <- rep(0, L)
    }
    
    new("immunoClust", expName=expName, parameters=parameters,
        K=L, P=P, N=N, w=obj$w[seq_len(L)], mu=mu, sigma=sigma,
        z=z, label=label, 
        logLike=obj$logLike, BIC=obj$logLike[1], ICL=obj$logLike[2],
        history=paste(obj$history), state=state)
    
}
## .immunoClust2