# Rajarshi Guha <rajarshi@presidency.com
# 5/12/2004
#
# Ref:Gao et al. JCICS, 1999, 39, 164-168
#     Labute, P., Pacific Symposium on Biocomputing 4:444-455 (1999)

require(KernSmooth)

.erf <- function(x) 2 * ( pnorm(x * sqrt(2)) ) -1

.get.params <- function(x,y, bintype='dpih') {
    # for each descriptor calculate bins
    # Uses dpih() from KernSmooth to get a bin width

    x1 <- x[ which(y==1), ]
    x0 <- x[ which(y==0), ]
    
    bwidth1 <- apply(x1, 2, dpih)
    bins1 <- list(); counts1 <- list()
    for (i in 1:ncol(x)) {
        b <- seq(min(x1[,i])-0.1, max(x1[,i])+0.1+bwidth1[i], by=bwidth1[i])
        c <- hist(x1[,i], breaks=b, plot=FALSE)$counts
        b[1] <- -Inf
        b[length(b)] <- Inf
        bins1[[i]] <- b
        counts1[[i]] <- c
    }
    
    bwidth0 <- apply(x0, 2, dpih)
    bins0 <- list(); counts0 <- list()
    for (i in 1:ncol(x)) {
        b <- seq(min(x0[,i])-0.1, max(x0[,i])+0.1+bwidth0[i], by=bwidth0[i])
        c <- hist(x0[,i], breaks=b, plot=FALSE)$counts
        b[1] <- -Inf
        b[length(b)] <- Inf
        bins0[[i]] <- b
        counts0[[i]] <- c
    }
    list(bins1=bins1, counts1=counts1, bins0=bins0, counts0=counts0)
}
   
.density <- function(z, bins, counts, sigma) {
    sum <- 0
    for (i in 2:length(bins)) {
        const <- sum(counts)
        tmp1 <- (counts[i-1] + 1/const)/(const + counts[i-1]/const)
        arg1 <- (bins[i] - z)/(sigma*sqrt(2))
        arg2 <- (bins[i-1] - z)/(sigma*sqrt(2))
        sum <- sum + tmp1 * (.erf(arg1) - .erf(arg2))
    }
    return(sum * 0.5)
}

# x data.frame of descriptors. It is a good idea to send in the rotated x
# onbtained by a PCA of the original data.
# y factor of class asignments (0,1), 1 indicating active
bqsar <- function(x,y, sigma=0.25) {

    m0 <- length(which(y == 0))
    m1 <- length(which(y == 1))
    print(paste("Pos = ",m1," Neg = ",m0,sep=""))

    params <- .get.params(x,y)
    attach(params)

    # for each descriptor calculate bins
    # Uses dpih() from KernSmooth to get a bind width
    p <- numeric(nrow(x))
    for (i in 1:nrow(x)) {
        prod <- 1
        for (j in 1:ncol(x)) {
            f0 <- .density(x[i,j], bins0[[j]], counts0[[j]], sigma)
            f1 <- .density(x[i,j], bins1[[j]], counts1[[j]], sigma)
            prod <- prod * f0/f1
        }
        p[i] <- 1.0 / (1.0 + (m0+1)*prod/(m1+1))
    }
    detach(params)
    p
}