# truecluster/R/information.R
# Some information theoretic functions
# (C) 2003-2007 Dr. Jens Oehlschlägel
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
###########################################################################

#source("c:/MWP/eAnalysis/truecluster/R/information.R")

#! \name{shannon.information}
#! \alias{laplace.probability}
#! \alias{laplace.information}
#! \alias{shannon.information}
#! \alias{absolute.entropy}
#! \alias{normalized.entropy}
#! \alias{redundancy}
#! \alias{distribution.entropy}
#! \alias{density.entropy}
#! \title{Univariate Informationtheoretic Measures}
#! \description{
#!   These functions calculate univariate informationtheoretic measures, like \emph{information}, \emph{entropy} and \emph{redundancy}.
#! }
#! \usage{
#! laplace.probability(n)
#! laplace.information(n)
#! shannon.information(p)
#! absolute.entropy(p)
#! normalized.entropy(p)
#! redundancy(p)
#! distribution.entropy(x, FUN = pnorm, ..., normalized = TRUE)
#! density.entropy(lower = -Inf, upper = Inf, FUN = dnorm, ..., normalized = TRUE)
#! }
#! \arguments{
#!   \item{n}{ scalar: number of (equal likely) possibilities }
#!   \item{p}{ vector of probabilities }
#!   \item{x}{ sequence of values at which FUN is evaluated }
#!   \item{lower}{ lower bound for integrate }
#!   \item{upper}{ upper bound for integrate }
#!   \item{FUN}{ a cumulative distribution function }
#!   \item{\dots}{ further arguments to FUN }
#!   \item{normalized}{ FALSE to use \code{absolute.entropy} instead of TRUE (default) \code{normalized.entropy} }
#! }
#! \value{
#!   Function \code{shannon.information} returns a vector of log2 (or zero where p=0). Entropy functions return vector of entropy contributions for each p.
#! }
#! \references{ MacKay, David J.C. (2003). \emph{Information Theory, Inference, and Learning Algorithms (chapter 8).} Cambridge University Press.}
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{twoway.entropy}}, \code{\link[distutil]{dist.entropy}}, \code{\link{Kullback.Leibler}}, \code{\link[base]{log}}, \code{\link{exact.margin.info}} }
#! \examples{
#!
#!   # elementary
#!   laplace.information(2)      # = log2(n)
#!   p <- laplace.probability(2) # = rep(1/n, n)
#!   p
#!   shannon.information(p)      # = log2(p)
#!   absolute.entropy(p)         # = p*log2(p)
#!   normalized.entropy(p)       # = p*log2(p) / log2(n)
#!
#!   sum(absolute.entropy(p))    # Max = log2(n)
#!   sum(normalized.entropy(p))  # Max = 1
#!   sum(redundancy(p))          # Redundancy = 1 - NormalizedEntropy
#!
#!   laplace.information(2)
#!   laplace.information(4)  # more categories, more information
#!   sum(absolute.entropy(laplace.probability(2)))
#!   sum(absolute.entropy(laplace.probability(4)))   # more categories, more entropy
#!   sum(normalized.entropy(laplace.probability(2)))
#!   sum(normalized.entropy(laplace.probability(4))) # more categories, constant normalized entropy
#!   p <- c(0.5, 0.25, 0.15, 0.1)
#!   sum(normalized.entropy(p))                      # unequal probabilities, lower entropy (normalized or not)
#!
#!   sum(distribution.entropy(seq(-3, 3, 0.01)))
#!   sum(distribution.entropy(seq(-3, 3, 0.001)))
#!   sum(distribution.entropy(seq(-30, 30, 0.01), sd=10))
#!   sum(distribution.entropy(seq(-3, 3, 0.01), FUN=punif, -3, 3))
#!
#!   sum(distribution.entropy(seq(-3, 3, 0.01), normalized=FALSE))
#!   sum(distribution.entropy(seq(-3, 3, 0.001), normalized=FALSE))
#!   sum(distribution.entropy(seq(-30, 30, 0.01), sd=10, normalized=FALSE))
#!   sum(distribution.entropy(seq(-3, 3, 0.01), FUN=punif, -3, 3, normalized=FALSE))
#!
#!   density.entropy(-3, 3)
#!   density.entropy(-30, 30, sd=10)
#!
#! }
#! \keyword{ univar }

laplace.probability <- function(n){
  rep(1/n, n)
}

laplace.information <- function(n){
  log(n, 2)
}

shannon.information <- function(p){
  i <- p>0
  I <- double(length(i))
  I[i] <- -log(p[i], 2)  # <-- this is the shannon formula, the rest is design decision
  I[is.na(I)] <- 0
  I
}

absolute.entropy <- function(p){
  p*shannon.information(p)
}

normalized.entropy <- function(p){
  H <- absolute.entropy(p)
  n <- sum(H>0)
  max.H <- laplace.information(n)
  H / max.H
}

redundancy <- function(p){
  H <- absolute.entropy(p)
  n <- length(p)
  max.H <- laplace.information(n)
  1/n - (H / max.H)
}

distribution.entropy <- function(x, FUN=pnorm, ..., normalized=TRUE){
  cp <- FUN(x, ...)
  p <- diff(cp)
  p <- p / sum(p)
  if (normalized){
    normalized.entropy(p)
  }else{
    warning("entropy increases with the number of categories")
    absolute.entropy(p)
  }
}

density.entropy <- function(lower=-Inf, upper=Inf, FUN=dnorm, ..., normalized=TRUE){
  warning("continuos entropy is meaningless: scale dependent")
  if (normalized){
    f <- function(x)normalized.entropy(FUN(x, ...))
    do.call("integrate", c(list(f=f, lower=lower, upper=upper)))
  }else{
    f <- function(x)absolute.entropy(FUN(x, ...))
    do.call("integrate", c(list(f=f, lower=lower, upper=upper)))
  }
}




#! \name{twoway.entropy}
#! \alias{twoway.entropy}
#! \alias{print.twoway.entropy}
#! \alias{plot.twoway.entropy}
#! \title{Bivariate Informationtheoretic Measures}
#! \description{
#!   These functions calculate bivariate informationtheoretic measures, like \emph{joint entropy}, \emph{conditional entropy} and \emph{mutual information}.
#! }
#! \usage{
#! twoway.entropy(p, grain=c("Totals","Margins","Cells"))
#! print.twoway.entropy(x, ...)
#! plot.twoway.entropy(x, ...)
#! }
#! \arguments{
#!   \item{p}{ matrix of two-way probabilities (or frequencies) }
#!   \item{x}{ on object of class twoway.entropy }
#!   \item{grain}{ one of c("Totals","Margins","Cells") determining the granularity of the return value }
#!   \item{\dots}{ further arguments to plot,print }
#! }
#! \value{
#!   An object of class twoway.entropy with components
#!   \item{Totals}{always returned: a list with components \emph{H} joint entropy, \emph{Ha} row entropy, \emph{Hab}  row conditional entropy given columns, \emph{Hb} column entropy, \emph{Hba} column conditional entropy given rows, \emph{Im} mutual information }
#!   \item{Margins}{returned unless grain="Total": \emph{pa} row probabilities , \emph{ha} row entropie, \emph{hab} colwise conditional entropies, \emph{pb} column probabilities, \emph{hb} column entropy , \emph{hba} rowwise conditional entropies}
#!   \item{Cells}{returned if grain="Cells": \emph{p} joint probabilities, \emph{pab} columnwise conditional probabilities, \emph{pba} rowwise conditional probabilities, \emph{h} joint entropies, \emph{hab} columnwise conditional entropies, \emph{hba} rowwise conditional entropies }
#! }
#! \references{ MacKay, David J.C. (2003). \emph{Information Theory, Inference, and Learning Algorithms (chapter 8).} Cambridge University Press.}
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{shannon.information}}, \code{\link[distutil]{dist.entropy}}, \code{\link{Kullback.Leibler}}, \code{\link[base]{log}} }
#! \examples{
#!  # Exercise 8.6 from MacKay
#!  p <- 1/matrix(c(8,16,16,4,16,8,16,Inf,32,32,16,Inf,32,32,16,Inf), 4)
#!  twoway.entropy(p)
#!  twoway.entropy(p, grain="Cells")
#!  str(twoway.entropy(p, grain="Cells"))
#!  plot(twoway.entropy(p), main="Plot from MacKay (2003), chapter 8.1")
#!
#!  # Highly correlated data
#!  x <- rnorm(100)
#!  y <- x + rnorm(100, sd=0.1)
#!  a <- round(x)
#!  b <- round(y)
#!  p <- table(a, b)
#!  plot(twoway.entropy(p), main=paste("correlation=", round(cor(x,y), 2)))
#! }
#! \keyword{ multivariate }


twoway.entropy <- function(p, grain=c("Totals","Margins","Cells")){
  grain <- match.arg(grain)
  if (!is.matrix(p))
    stop("x must be matrix")
  p <- p / sum(p)
  pa <- rowSums(p)
  pb <- colSums(p)

  h  <- absolute.entropy(p)
  ha <- absolute.entropy(pa)
  hb <- absolute.entropy(pb)

  H <- sum(h)
  Ha <- sum(ha)
  Hb <- sum(hb)
  Hab <- H - Hb
  Hba <- H - Ha
  Im <- H - Hab - Hba

  if (grain!="Totals"){
    pab <- t(t(p)/pb)
    pba <- p/pa
    hab. <- apply(pab, 2, absolute.entropy)
    hba. <- t(apply(pba, 1, absolute.entropy))
    hab <- colSums(hab.)
    hba <- rowSums(hba.)
  }

  ret <- list()
  if (grain=="Cells"){
    ret[["Cells"]] <- list(p=p, pab=pab, pba=pba, h=h, hab=hab., hba=hba.)
  }
  if (grain!="Totals"){
    ret[["Margins"]] <- list(pa=pa, ha=ha, hab=hab, pb=pb, hb=hb, hba=hba)
  }
  ret[["Totals"]] <- list(H=H, Ha=Ha, Hab=Hab, Hb=Hb, Hba=Hba, Im=Im)
  class(ret) <- "twoway.entropy"
  ret
}

print.twoway.entropy <- function(x, ...){
  Totals <- unlist(x$Totals)
  if (is.null(x$Margins)){
    Margins <- NULL
  }else{
    m <- x$Margins
    Ha  <- sum(m$ha)
    Hb  <- sum(m$hb)
    Hab <- sum(m$pb*m$hab)
    Hba <- sum(m$pa*m$hba)
    Im  <- Ha-Hab
    H   <- Hab + Hba + Im
    Margins <- c(H=H, Ha=Ha, Hab=Hab, Hb=Hb, Hba=Hba, Im=Im)
  }
  print(rbind("via Totals"=Totals, "via Margins"=Margins), ...)
}


plot.twoway.entropy <- function(x, ...){
  m <- x$Totals
  dat <- cbind(
    Mutual  = c(Hab=m$Hab, Im=m$Im, Hba=m$Hba)
  , Columns = c(Hab=m$Hab, " "=0, Hb=m$Hb)
  , Rows    = c(Ha=m$Ha, " "=0, Hba=m$Hba)
  , Total   = c(H=m$H, " "=0, " "=0)
  )
  mid <- barplot(dat, horiz=TRUE, col="transparent", ...)
  dat2 <- apply(rbind(rep(0,4), dat), 2, cumsum)
  dat3 <- (dat2[-1,] + dat2[-4,])/2
  text(dat3[,"Mutual"], mid[1], c("Hab","Im","Hba"))
  text(dat3[,"Columns"], mid[2], c("Hab","","Hb"))
  text(dat3[,"Rows"], mid[3], c("Ha","","Hba"))
  text(dat3[,"Total"], mid[4], c("H","",""))
}





#! \name{dist.entropy}
#! \alias{dist.entropy}
#! \title{ Entropy distance between variables }
#! \description{
#!   This functions calculates the \emph{entropy distance} between variables
#! }
#! \usage{
#! dist.entropy(x)
#! }
#! \arguments{
#!   \item{x}{ a matrix of categorical raw data (cases in rows, variables in columns) }
#! }
#! \value{
#!   An object of class \code{\link[stat]{dist}}.
#! }
#! \references{
#!   MacKay, David J.C. (2003). \emph{Information Theory, Inference, and Learning Algorithms (Exercise 8.5).} Cambridge University Press.
#! }
#! \author{ Jens Oehlschlaegel }
#! \seealso{ \code{\link[stat]{dist}}, \code{\link[rjens]{cor2dist}}, \code{\link[rjens]{cor2sim}}, \code{\link[rjens]{sim2dist}} }
#! \examples{
#!   n <- 1000
#!   x <- runif(n)
#!   x <- cbind(
#!     x
#!   , x
#!   , 2*x
#!   , x+runif(n, max=0.1)
#!   , x+runif(n)
#!   , x+runif(n, max=10)
#!   , x+runif(n, max=10)
#!   , 2*(x+runif(n, max=10))
#!   , 2*(x+runif(n, max=10))
#!   , runif(n, max=22)
#!   , runif(n, max=22)
#!   )
#!   x <- round(x)
#!   colnames(x) <- letters[1:ncol(x)]
#!   D <- dist.entropy(x)
#!   plot(hclust(D, method="average"))
#! }
#! \keyword{multivariate}
#! \keyword{cluster}

dist.entropy <- function(x){
  n <- ncol(x)
  nam <- colnames(x)
  D <- matrix(as.double(NA), n, n, dimnames=list(nam, nam))
  for (i in 2:n)
  for (j in 1:(i-1)){
    te <- twoway.entropy(as.matrix(table(x[,c(i,j)])))
    D[i,j] <- D[j,i] <- te$Totals$H - te$Totals$Im
  }
  as.dist(D)
}


#! \name{Kullback.Leibler}
#! \alias{Kullback.Leibler}
#! \title{Multivariate Informationtheoretic Measures}
#! \description{
#!   These functions calculate mulivariate informationtheoretic measures: \emph{Kullback Leibler divergence}.
#! }
#! \usage{
#! Kullback.Leibler(p, basis = 2)
#! }
#! \arguments{
#!   \item{p}{ matrix, each column contains probabilities of one distribution across the same alphabet (rows) }
#!   \item{basis}{ basis of the log used (default=2) }
#! }
#! \details{
#!   \deqn{D_{KL}(P||Q) = \sum(P(x) \dot log \frac{P(x)}{Q(x)})}{D_KL(P||Q) = sum(P*log(P/Q))}
#! }
#! \value{
#!   A square divergence matrix,
#!   \item{Totals}{always returned: a list with components \emph{H} joint entropy, \emph{Ha} row entropy, \emph{Hab}  row conditional entropy given columns, \emph{Hb} column entropy, \emph{Hba} column conditional entropy given rows, \emph{Im} mutual information }
#!   \item{Margins}{returned unless grain="Total": \emph{pa} row probabilities , \emph{ha} row entropie, \emph{hab} colwise conditional entropies, \emph{pb} column probabilities, \emph{hb} column entropy , \emph{hba} rowwise conditional entropies}
#!   \item{Cells}{returned if grain="Cells": \emph{p} joint probabilities, \emph{pab} columnwise conditional probabilities, \emph{pba} rowwise conditional probabilities, \emph{h} joint entropies, \emph{hab} columnwise conditional entropies, \emph{hba} rowwise conditional entropies }
#! }
#! \references{ MacKay, David J.C. (2003). \emph{Information Theory, Inference, and Learning Algorithms (chapter 8).} Cambridge University Press.}
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{shannon.information}}, \code{\link[distutil]{dist.entropy}}, \code{\link{Kullback.Leibler}}, \code{\link[base]{log}} }
#! \examples{
#!  x <- seq(-3, 3, 0.1)
#!  cp <- pnorm(x)
#!  p <- cp[-1] - cp[-length(cp)]
#!  cq <- punif(x, -3, 3)
#!  q <- cq[-1] - cq[-length(cq)]
#!  Kullback.Leibler(cbind(p,q))
#! }
#! \keyword{ multivariate }

Kullback.Leibler <- function(p, basis=2){
  if (!is.matrix(p))
    stop("need matrix p of column probabilities")
  margin <- colSums(p)
  if (any(margin!=1)){
    warning("probabilitites do not sum to 1 colwise, will normalize")
    p <- t(t(p) / margin)
  }
  if (any(p<=0))
    stop("not all p positive")
  l <- log(p,basis)
  n <- ncol(p)
  nam <- colnames(p)
  D <- matrix(as.double(0), n, n, dimnames=list(nam,nam))
  for (i in 1:n)
  for (j in 1:i){
    dl <- l[,i]-l[,j]
    D[i,j] <- sum(p[,i]*dl)
    D[j,i] <- -sum(p[,j]*dl)
  }
  D
}

