# truecluster/R/mmcc.R
# Truecluster CIC evaluation of MMCC results
# (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/mmcc.r")

require(distutil)  # because of rowMax (since max.col does certain cases wrong in R versions prior 2.5.0)

#! \name{mmcc2mmcp}
#! \alias{mmcc2mmcp}
#! \alias{mmcc2cluster}
#! \alias{mmcp2pk}
#! \alias{mmcp2RMC}
#! \alias{mmcp2h}
#! \alias{mmcp2H}
#! \alias{mmcp2d}
#! \alias{mmcp2D}
#! \alias{mmcp2i}
#! \alias{mmcp2I}
#! \alias{mmcp2cic}
#! \alias{mmcp2CIC}
#! \alias{mmcp2ALL}
#! \alias{mmcp2gsd}
#! \alias{mmcp2GSD}
#! \alias{cic}
#! \alias{CIC}
#! \alias{mmcc}
#! \alias{MMCC}
#! \alias{gsd}
#! \alias{GSD}
#! \title{MMCC Evaluation}
#! \description{
#!   These functions evaluate the MMCC matrix created by resample aggregation.
#! }
#! \usage{
#! mmcc2cluster(c)
#! mmcc2mmcp(c, rowsums = NULL, equal.rowsums = FALSE)
#! mmcp2pk(p)   # marginal cluster probabilities
#! mmcp2RMC(p)  # relative model complexity
#! mmcp2h(p)    # cellwise conditional entropy = uncertainty
#! mmcp2H(p)    # conditional entropy = uncertainty
#! mmcp2d(p)    # cellwise non-penalized information
#! mmcp2D(p)    # non-penalized information
#! mmcp2i(p)    # cellwise penalized information
#! mmcp2I(p)    # penalized information
#! mmcp2cic(p)  # cellwise CIC
#! mmcp2CIC(p)  # CIC
#! mmcp2ALL(p)  # c(N, H, D, RMC, I, CIC)
#! mmcp2gsd(p)  # rowwise Generalized Silhouette Value
#! mmcp2GSD(p, FUN=summary)  # clusterwise statistics of Generalized Silhouette Value
#! }
#! \arguments{
#!   \item{c}{MMCC matrix of cluster votes}
#!   \item{p}{MMCP matrix of cluster probabilities}
#!   \item{rowsums}{specifying vector or scalar of rowsums of mmcc speeds up calculations}
#!   \item{equal.rowsums}{TRUE uses sum(mmcc[1,]) as rowsums, ignored if rowsums is given }
#!   \item{FUN}{a function summarizing a vector of GSD values (called once per cluster)}
#! }
#! \details{
#!  \code{mmcc2cluster} assigns each case (row) to the most likely cluster (column) breaking ties at random.
#! }
#! \value{
#!   Function mmcc2cluster returns a vector of cluster memberships.
#!   Function mmcc2mmcp returns converts an MMCC (counted cluster votes) matrix into an MMCP matrix (estimated cluster probabilities).
#!   Function mmcp2RMC returns scalar giving the \emph{relative model complexity}
#!   Function mmcp2gsd returns a vector of Generalized Silhouette Diagnostics
#!   Function mmcp2GSD returns a matrix[statistics,clusters] summarizing Generalized Silhouette Diagnostics
#!   Function mmcp2xx returns a matrix derived from MMCP.
#!   Function mmcp2XX returns mean(rowSums(mmcp2xx())) and might be a bit faster.
#! }
#! \references{ xx }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{twoway.entropy}}, \code{\link{ldmargination}} }
#! \examples{
#! cat("\nfake a cluster voting matrix\n")
#! c <- c(rbinom(4, 256, 0.25), rbinom(12, 256, 0.75))
#! c <- cbind(a=c, b=256-c)
#! c
#! mmcc2cluster(c)
#! table(mmcc2cluster(c))
#! cat("\nconvert to probabilities matrix\n")
#! p <- mmcc2mmcp(c)
#! p
#!
#! cat("\nlook at Relative Model Complexity\n")
#! mmcp2RMC(p)
#!
#! cat("\ncellwise Cluster Information Criterion\n")
#! mmcp2cic(p)
#! cat("\ncluster CIC\n")
#! colSums(mmcp2cic(p))
#! cat("\ncase CIC (aggregated by cluster)\n")
#! aggregate(rowSums(mmcp2cic(p)), by=list(mmcc2cluster(p)), FUN=mean)
#! cat("\ncalculating model CIC\n")
#! mean(rowSums(mmcp2cic(p)))
#! mmcp2CIC(p)
#! cat("\ncalculating all model statistics\n")
#! mmcp2ALL(p)
#! cat("\nGeneralized Silhouette Diagnostic\n")
#! mmcp2GSD(p)
#! boxplot(split(mmcp2gsd(p), mmcc2cluster(p)), main="GSD distribution by cluster")
#! boxplot(split(rowSums(mmcp2h(p)), mmcc2cluster(p)), main="Uncertainty distribution by cluster")
#!
#! cat("\nshow that ties are broken at random\n")
#! c <- matrix(1, 16, 4)
#! table(mmcc2cluster(c))
#! table(mmcc2cluster(c))
#!
#! cat("\ncompare some models differing in true and model complexity\n")
#! collect <- NULL
#! c <- matrix(1, 16, 1)
#! collect <- rbind(collect, "true=1  model=1"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=1  model=2"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=1  model=4"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=1  model=8"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=1  model=16"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- t(matrix(rep(diag(2), 8), 2))
#! collect <- rbind(collect, "true=2  model=2"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=2  model=4"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=2  model=8"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=2  model=16"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- t(matrix(rep(diag(4), 4), 4))
#! collect <- rbind(collect, "true=4  model=4"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=4  model=8"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=4  model=16"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- t(matrix(rep(diag(8), 2), 8))
#! collect <- rbind(collect, "true=8  model=8"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- cbind(c, c)
#! collect <- rbind(collect, "true=8  model=16"=mmcp2ALL(mmcc2mmcp(c)))
#! c <- diag(16)
#! collect <- rbind(collect, "true=16 model=16"=mmcp2ALL(mmcc2mmcp(c)))
#! colnames(collect)[c(5,6)] <- c("D*(1-RMC)=I","I-H=CIC")
#! print(collect[order(collect[,6]),])
#!
#! }
#! \keyword{multivariate}
#! \keyword{cluster}


# xx could use rowMax directly instead of
mmcc2cluster <- function(c)rowMax(c, value=FALSE)

mmcc2mmcp <- function(
  c
, rowsums          = NULL     # scalar or vector to speed up calculations
, equal.rowsums    = FALSE    # TRUE to calculate rowsums just as scalar (ignored if rowsums given)
)
  # RETURNs matrix[cases, clusters]
{
  if (is.null(rowsums)){
    if (equal.rowsums){
      rowsums <- sum(c[1, ])
    }else{
      rowsums <- rowSums(c)
    }
  }
  c / rowsums
}


mmcp2pk <- function(p){
  colMeans(p)
}

mmcp2RMC <- function(p){
  pk <- colMeans(p)
  hk <- -pk*log(pk, 2)
  hk[pk==0] <- 0
  (2^sum(hk) - 1) / (nrow(p)-1)
}

mmcp2h <- function(p){
  h <- -p*log(p, 2)
  h[p==0] <- 0
  h
}

mmcp2H <- function(p){
  h <- -p*log(p, 2)
  h[p==0] <- 0
  sum(h)/nrow(p)
}

mmcp2d <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    p[] <- 0
    p
  }else{
    -p*log(1-abs(t(t(p)-pk)), 2)
  }
}

mmcp2D <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    0
  }else{
    sum(-p*log(1-abs(t(t(p)-pk)), 2))/nrow(p)
  }
}

mmcp2i <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    p[] <- 0
    p
  }else{
    hk <- -pk*log(pk, 2)
    hk[pk==0] <- 0
    relative.complexity <- (2^sum(hk) - 1) / (nrow(p)-1)
    D <- -p*log(1-abs(t(t(p)-pk)), 2)
    D * (1-relative.complexity)
  }
}

mmcp2I <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    0
  }else{
    n <- nrow(p)
    hk <- -pk*log(pk, 2)
    hk[pk==0] <- 0
    relative.complexity <- (2^sum(hk) - 1) / (n-1)
    D <- sum(-p*log(1-abs(t(t(p)-pk)), 2)) / n
    D * (1-relative.complexity)
  }
}



mmcp2cic <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    p[] <- 0
    p
  }else{
    hk <- -pk*log(pk, 2)
    hk[pk==0] <- 0
    relative.complexity <- (2^sum(hk) - 1) / (nrow(p)-1)
    D <- -p*log(1-abs(t(t(p)-pk)), 2)
    I <- D * (1-relative.complexity)
    H <- p
    i <- as.vector(p>0)
    H[i] <- -p[i]*log(p[i], 2)
    I - H
  }
}

mmcp2CIC <- function(p){
  pk <- colMeans(p)
  if (any(pk==1)){
    0
  }else{
    n <- nrow(p)
    hk <- -pk*log(pk, 2)
    hk[pk==0] <- 0
    relative.complexity <- (2^sum(hk) - 1) / (n-1)
    D <- sum(-p*log(1-abs(t(t(p)-pk)), 2)) # /n deferred to later
    I <- D * (1-relative.complexity)
    H <- p
    i <- as.vector(p>0)
    H[i] <- -p[i]*log(p[i], 2)
    (I - sum(H))/n
  }
}

mmcp2ALL <- function(p){
  n <- nrow(p)
  pk <- colMeans(p)
  hk <- -pk*log(pk, 2)
  hk[pk==0] <- 0
  relative.complexity <- (2^sum(hk) - 1) / (nrow(p)-1)
  D <- -sum(p*log(1-abs(t(t(p)-pk)), 2))
  I <- D * (1-relative.complexity)
  if (any(pk==1)){
    D <- 0
    I <- 0
  }
  H <- p
  i <- as.vector(p>0)
  H[i] <- -p[i]*log(p[i], 2)
  H <- sum(H)
  c(N=n, H=H/n, D=D/n, RMC=relative.complexity, I=I/n, CIC=(I-H)/n)
}

mmcp2gsd <- function(p){
  n <- nrow(p)
  i <- cbind(seq(length=n), rowMax(p, value=FALSE))
  p1 <- p[i]
  p[i] <- -Inf
  p2 <- rowMax(p)
  (p1/(p1+p2)) * 2 - 1
}

mmcp2GSD <- function(p, FUN=summary){
  n <- nrow(p)
  i <- cbind(seq(length=n), rowMax(p, value=FALSE))
  p1 <- p[i]
  p[i] <- -Inf
  p2 <- rowMax(p)
  gsd <- (p1/(p1+p2)) * 2 - 1
  sapply(split(gsd, i[,2]), FUN)
}
