# truecluster/R/combinatorics.R
# Some combinatorial 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/combinatorics.r")

#! \name{combination}
#! \alias{combination}
#! \alias{ncombination}
#! \title{ Combinations }
#! \description{
#!   Calculating and creating combinations
#! }
#! \usage{
#! ncombination(k, n, replace = FALSE)
#! combination(k, n=NULL, x=NULL, replace = FALSE)
#! }
#! \arguments{
#!   \item{k}{ number of positions }
#!   \item{n}{ poolsize }
#!   \item{x}{ vector pool of values, can be given instead of n }
#!   \item{replace}{ default FALSE, TRUE to allow reusing of pool elements }
#! }
#! \details{
#!   \code{ncombination} returns the number of combinations, \code{combination} returns a matrix of combinations. No special sequence granted. Implemented in C, but not faster than the functions in gregmisc.
#! }
#! \value{
#!   A matrix of (k x ncombination)
#! }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{permutation}}, \code{\link{labelation}}, \code{\link[base]{choose}}, \code{\link[gregmisc]{permutations}} }
#! \examples{
#!   ncombination(2,4)
#!   combination(2,4)
#!   combination(2,x=letters[1:4])
#!   ncombination(2,4,replace=TRUE)
#!   combination(2,4,replace=TRUE)
#!
#!   \dontrun{
#!   ## regression test
#!   require(gregmisc)
#!   k <- 4
#!   n <- 6
#!   replace <- FALSE
#!   a <- sort(apply(combination(k,n,replace=replace), 2, paste, collapse="-"))
#!   b <- sort(apply(combinations(n,k,rep=replace), 1, paste, collapse="-"))
#!   stopifnot(length(a)==ncombination(k,n,replace=replace))
#!   stopifnot(identical(a,b))
#!   replace <- TRUE
#!   a <- sort(apply(combination(k,n,replace=replace), 2, paste, collapse="-"))
#!   b <- sort(apply(combinations(n,k,rep=replace), 1, paste, collapse="-"))
#!   stopifnot(length(a)==ncombination(k,n,replace=replace))
#!   stopifnot(identical(a,b))
#!   }
#! }
#! \keyword{ iteration }

ncombination <- function(k, n, replace=FALSE){
  if (replace){
    x <- 1:n
    for (i in 2:k)
      x <- cumsum(x)
    x[n]
  }else{
    choose(n, k)
  }
}

combination <- function (k, n=NULL, x=NULL, replace=FALSE)
{
    if (is.null(n)){
      n <- length(x)
      if (!length(n))
        stop("need length(x)>0 when is.null(n)")
    }else{
      if (!is.null(x) && length(x)!=n)
        stop("if n and x is given we need n=length(x)")
    }
    if (k>n && !replace)
      stop("k>n not possible with replace=FALSE")
    C <- ncombination(k, n, replace=replace)
    if (is.na(as.integer(C*k)))
      stop("too big", C)
    if (replace){
      ret <- .C("R_tc_combination_replace", M=integer(C*k), N=as.integer(n), K=as.integer(k),
          NAOK = FALSE, DUP = FALSE, PACKAGE = "truecluster")
    }else{
      ret <- .C("R_tc_combination_noreplace", M=integer(C*k), N=as.integer(n), K=as.integer(k),
          NAOK = FALSE, DUP = FALSE, PACKAGE = "truecluster")
    }
    if (is.null(x))
      matrix(ret$M+1, nrow=k)
    else
      matrix(x[ret$M+1], nrow=k)
}


#! \name{permutation}
#! \alias{permutation}
#! \alias{npermutation}
#! \title{ Permutations }
#! \description{
#!   Calculating and creating permutations
#! }
#! \usage{
#! npermutation(k, n, replace = FALSE)
#! permutation(k, n=NULL, x=NULL, replace = FALSE)
#! }
#! \arguments{
#!   \item{k}{ number of positions }
#!   \item{n}{ poolsize }
#!   \item{x}{ vector pool of values, can be given instead of n }
#!   \item{replace}{ default FALSE, TRUE to allow reusing of pool elements }
#! }
#! \details{
#!   \code{npermutation} returns the number of permutations, \code{permutation} returns a matrix of permutations. No special sequence granted. Implemented in C, but not faster than the functions in gregmisc.
#! }
#! \value{
#!   A matrix of (k x npermutation)
#! }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{combination}}, \code{\link{margination}}, \code{\link[base]{expand.grid}}, \code{\link[gregmisc]{combinations}} }
#! \examples{
#!   npermutation(2,4)
#!   permutation(2,4)
#!   permutation(2,x=letters[1:4])
#!   npermutation(2,4,replace=TRUE)
#!   permutation(2,4,replace=TRUE)
#!
#!   \dontrun{
#!   ## regression test (error in gregmisc for 1.5.1)
#!   require(gregmisc)
#!   k <- 3
#!   n <- 6
#!   replace <- FALSE
#!   a <- sort(apply(permutation(k,n,replace=replace), 2, paste, collapse="-"))
#!   b <- sort(apply(permutations(n,k,rep=replace), 1, paste, collapse="-"))
#!   stopifnot(length(a)==npermutation(k,n,replace=replace))
#!   stopifnot(identical(a,b))
#!   replace <- TRUE
#!   a <- sort(apply(permutation(k,n,replace=replace), 2, paste, collapse="-"))
#!   b <- sort(apply(permutations(n,k,rep=replace), 1, paste, collapse="-"))
#!   stopifnot(length(a)==npermutation(k,n,replace=replace))
#!   stopifnot(identical(a,b))
#!   }
#! }
#! \keyword{ iteration }

npermutation <- function(k, n, replace=FALSE){
  if (replace){
    n^k
  }else{
    eval(parse(text=paste(seq(n, by=-1, length=k), collapse="*")))
  }
}


permutation <- function (k, n=NULL, x=NULL, replace=FALSE)
{
    if (is.null(n)){
      n <- length(x)
      if (!length(n))
        stop("need length(x)>0 when is.null(n)")
    }else{
      if (!is.null(x) && length(x)!=n)
        stop("if n and x is given we need n=length(x)")
    }
    if (k>n && !replace)
      stop("k>n not possible with replace=FALSE")
    C <- npermutation(k, n, replace=replace)
    if (is.na(as.integer(C*k)))
      stop("too big", C)
    if (replace){
      ret <- .C("R_tc_permutation_replace", M=integer(C*k), N=as.integer(n), K=as.integer(k),
          NAOK = FALSE, DUP = FALSE, PACKAGE = "truecluster")
    }else{
      ret <- .C("R_tc_permutation_noreplace", M=integer(C*k), N=as.integer(n), K=as.integer(k),
          NAOK = FALSE, DUP = FALSE, PACKAGE = "truecluster")
    }
    if (is.null(x))
      matrix(ret$M+1, nrow=k)
    else
      matrix(x[ret$M+1], nrow=k)
}


#! \name{margination}
#! \alias{labelation}
#! \alias{nlabelation}
#! \alias{lnlabelation}
#! \alias{ldlabelation}
#! \alias{margination}
#! \alias{nmargination}
#! \alias{lnmargination}
#! \alias{ldmargination}
#! \title{ Margin combinations }
#! \description{
#!   Functions to calculate the number of possible groupings given fixed group sizes (margins).
#! }
#! \usage{
#! nlabelation(...)
#! lnlabelation(...)
#! ldlabelation(...)
#! nmargination(...)
#! lnmargination(...)
#! ldmargination(...)
#! }
#! \arguments{
#!   \item{\dots}{ one or more vectors of group sizes }
#! }
#! \details{
#!   Function \code{nmargination} gives the number of possibles groupings (=cluster solutions) that have the input group sizes.
#!   Function \code{nlabelation} gives the number of possibles assignents to predefined (=labelled) groups of the requires sizes (more possibilities).
#!   The ln-versions give the log instead, the ld-versions the logarithm dualis.
#! }
#! \value{
#!   scalar
#! }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{clusteration}}, \code{\link{permutation}}, \code{\link{combination}}, \code{\link[base]{choose}}, \code{\link[base]{lgamma}}, \code{\link{absolute.entropy}}, \code{\link{exact.margin.info}} }
#! \examples{
#!   cat("There are 25 possibilities to make 3 clusters from 5 cases\n")
#!   nclusteration(5,3)
#!   cat(" 10 of sizes 3,1,1\n")
#!   nmargination(3,1,1)
#!   cat(" and 15 of sizes 2,2,1\n")
#!   nmargination(2,2,1)
#!   cat("By contrast there are 30 possibilities to assign 5 cases to 3 labelled groups a=2,b=2,c=1\n")
#!   nlabelation(2,2,1)
#!   n <- 250
#!   k <- 5
#!   i <- 1:(n/k)
#!   matplot(k*i, cbind(
#!       k^(i*k)
#!     , sapply(i, function(i){nclusteration(k*i, k)})
#!     , sapply(i, function(i){nmargination(rep(i, k))})
#!     )
#!   , log="y", type="l"
#!   , lty=1:3, xlab="n"
#!   , ylab = "combinations"
#!   , main = "number of possible groupings"
#!   )
#!   legend(n,1
#!   , c(substitute(k^n, list(k=k)), expression("no restriction"), expression("equal sizes"))
#!   , lty=1:3, col=1:3, xjust=1, yjust=0
#!   )
#! }
#! \keyword{ math }
#! \keyword{ cluster }


nlabelation <- function(...)
{
  tab <- c(...)
  k <- length(tab)
  if (k>2){
    i <- k%/%2
    choose(sum(tab), sum(tab[1:i])) * nlabelation(tab[1:i]) * nlabelation(tab[(i+1):k])
  }else if (k==2){
    choose(sum(tab), tab[1])
  }else 1
}

lnlabelation <- function(...)
{
  tab <- c(...)
  lgamma(sum(tab)+1) - sum(lgamma(tab+1))
}
ldlabelation <- function(...){
  tab <- c(...)
  lnlabelation(tab) / log(2)
}


nmargination <- function(...)
{
  tab <- c(...)

  # determine chunks of clusters with same size
  chunks <- table(tab)
  chunks <- cbind(as.double(names(chunks)), as.double(chunks), as.double(names(chunks)) * as.double(chunks))

  # identify singletons
  # and find out no. of possibilities to put cases into singleton set
  chunks <- chunks[order(chunks[,1]), , drop=FALSE]
  if (chunks[1,1]==1){
    singleposs <- choose(sum(chunks[,3]), chunks[1,3])
    chunks <- chunks[-1, , drop=FALSE]
  }else{
    singleposs <- 1
  }
  nchunks <- nrow(chunks)
  if (nchunks){
    # find out no. of possibilities to put non-singleton cases into chunks
    # as chunks have their meaning externally by size we can use nlabelation
    chunkposs <- nlabelation(chunks[,3])
    # find out no. of possibilities within each chunk
    innerposs <- sapply(1:nchunks, function(i){
      tab <- rep(chunks[i, 1], chunks[i, 2])
      prod(choose(cumsum(tab)-1, tab-1))
    })
  }else{
    chunkposs <- 1
    innerposs <- 1
  }
  prod(c(singleposs, chunkposs, innerposs))
}


lnmargination <- function(...)
{
  tab <- c(...)

  # determine chunks of clusters with same size
  chunks <- table(tab)
  if (any(names(chunks)=="0")){
    #warning("removing groups of size zero")
    chunks <- chunks[names(chunks)!="0"]
  }
  sizes <- as.integer(names(chunks))
  freqs <- as.integer(chunks)
  # find out no. of possibilities to split cases into chunks
  # as chunks have their meaning externally by size we can use lnlabelation
  lchunkposs <- lnlabelation(sizes*freqs)
  # find out no. of possibilities within each chunk
  linnerposs <- sapply(1:length(chunks), function(i){
    tab <- rep(sizes[i], freqs[i])
    sum(lchoose(cumsum(tab)-1, tab-1))
  })
  sum(c(lchunkposs, linnerposs))
}

ldmargination <- function(...){
  lnmargination(...) / log(2)
}


#! \name{clusteration}
#! \alias{clusteration}
#! \alias{nclusteration}
#! \alias{print.clusteration}
#! \title{ Cluster combinations }
#! \description{
#!   Functions to construct all possible groupings or calculate the number of possible groupings.
#! }
#! \usage{
#! clusteration(n = NULL, k = NULL, nmax = NULL, kmax = nmax)
#! nclusteration(n = NULL, k = NULL, nmax = NULL, kmax = nmax)
#! \method{print}{clusteration}(x, within = "-", between = "+", \dots)
#! }
#! \arguments{
#!   \item{n}{ number of elements }
#!   \item{k}{ number of groups (optional) }
#!   \item{nmax}{ maximal number of elements }
#!   \item{kmax}{ maximal number of groups }
#!   \item{x}{ an object of class 'clusteration' }
#!   \item{within}{ character to seperate elements within the same grouping }
#!   \item{between}{ charachter to seperate between groupings }
#!   \item{\dots}{ ignored }
#! }
#! \details{
#!   Can be used to create/calculate all groupings (using \code{nmax}), where the number of groupings can be restricted (\code{kmax}).
#!   The result is a matrix (n x k) which contains in the \code{\link[base]{lower.tri}} Jens' triangle not unlike Pascal's triangle.
#!   \cr
#!   Using \code{n} instead of \code{nmax} will return a vector equaling the last line of the triangle,
#!   using \code{k} can be used to directly return the groupings/number of groupings for a single fixed number of groups.
#!   All results are created using the triangle construction method
#!   just calling \code{nclusteration(n=, k=) makes use of Jensen's formula}
#!   \deqn{C(n,k)=\frac{1}{k!} \cdot \sum_{c=1}^{k}{\over{k}{c} \cdot k^n }}
#!   The triangle construction method uses the recursive formula
#!   \deqn{C(n,k)=C(n-1,k-1) + k \cdot C(n-1,k)}
#!   where \eqn{C(n,1)=1} and \eqn{C(n,n)=1}.
#! }
#! \value{
#!   Function \code{nclusteration} returns either a matrix[n,k] or a vector[k] or a scalar.
#!   Function \code{clusteration} returns either a list matrix[[n,k]] or a list vector[[k]] or a single element list vector[[1]] of list of \emph{groupings},
#!   where \emph{groupings} is a list with as many elements as groups and each group beeing represented by a vector of group member indices. A special print method prints this nicely.
#! }
#! \references{ Jensen, R.E. (1969) A dynamic programming algorithm for cluster analysis. \emph{Operations Research}, 17, 1034-57. }
#! \author{ Jens Oehlschlägel }
#! \note{ The total number of possible groupings increases worse than exponential }
#! \seealso{ \code{\link{margination}}, \code{\link{combination}}, \code{\link[base]{choose}} }
#! \examples{
#!    clusteration(nmax=5)
#!    clusteration(n=5)
#!    clusteration(n=5, k=3)
#!    nclusteration(nmax=5)
#!    nclusteration(n=5)
#!    nclusteration(n=5, k=3)
#!
#!    n <- 100
#!    k <- 5
#!    x <- nclusteration(nmax=n)
#!    matplot(1:n, cbind((1:n)^(1:n), rowSums(x), k^(1:n), x[,k])
#!    , log="y", type="l"
#!    , lty=1:4, xlab="n"
#!    , ylab = "combinations"
#!    , main = "number of possible groupings"
#!    )
#!    legend(n,1
#!    , c(expression(n^n), expression(k \%in\% group("{", list(1,2,...,n), "}")), substitute(k^n, list(k=k)), substitute(k==c, list(c=k))
#!    ), lty=1:4, col=1:4, xjust=1, yjust=0)
#! }
#! \keyword{ math }
#! \keyword{ cluster }


nclusteration <- function(n=NULL,k=NULL, nmax=NULL, kmax=nmax){
  if (is.null(nmax)){
    stopifnot(!is.null(n))
    if (is.null(k)){
      if (is.null(kmax))
        kmax <- n
      m <- rep(as.integer(0), kmax)
      names(m) <- 1:kmax
      for (n in 1:n){
        pre <- m
        for (k in 1:min(n,kmax)){
          if (k==1 || k==n)
            m[k] <- 1
          else
            m[k] <- pre[k-1] + k*pre[k]
        }
      }
    }else{
      # Jensen 1969, according to Gordon 1999, 2nd edition, p. 40
      c <- 1:k
      m <- sum((-1)^(k-c)*choose(k,c)*c^n)/(gamma(k+1)) # factorial(k) not available in all versions
    }
  }else{
    stopifnot(is.null(n) && is.null(k))
    m <- matrix(as.integer(0), nrow=nmax, ncol=kmax, dimnames=list(n=1:nmax, k=1:kmax))
    for (n in 1:nmax){
      for (k in 1:min(n,kmax)){
        if (k==1 || k==n)
          m[n,k] <- 1
        else
          m[n,k] <- m[n-1,k-1] + k*m[n-1,k]
      }
    }
  }
  m
}


clusteration <- function(n=NULL,k=NULL, nmax=NULL, kmax=nmax){
  if (is.null(nmax)){
    stopifnot(!is.null(n))
    if (is.null(kmax))
      kmax <- n
    m <- vector(mode="list", kmax)
    names(m) <- 1:kmax
    kk <- k
    for (n in 1:n){
      pre <- m
      for (k in 1:min(n,kmax)){
        if (k==1)
          m[[k]] <- list(list(1:n))
        else if (k==n)
          m[[k]] <- list(as.list(1:n))
        else{
          m[[k]] <- c(
            lapply(pre[[k-1]], c, n)
          , unlist(lapply(pre[[k]], function(l){
              lapply(1:length(l), function(i){
                l2 <- l
                l2[[i]] <- c(l2[[i]], n)
                l2
              })
            }), recursive=FALSE)
          )
        }
      }
    }
    if (!is.null(kk))
      m <- m[kk]
    attr(m, "n") <- n
  }else{
    stopifnot(is.null(n) && is.null(k))
    m <- matrix(list(NULL), nrow=nmax, ncol=kmax, dimnames=list(n=1:nmax, k=1:kmax))
    for (n in 1:nmax){
      for (k in 1:min(n,kmax)){
        if (k==1)
          m[[n,k]] <- list(list(1:n))
        else if (k==n)
          m[[n,k]] <- list(as.list(1:n))
        else{
          m[[n,k]] <- c(
            lapply(m[[n-1,k-1]], c, n)
          , unlist(lapply(m[[n-1,k]], function(l){
              lapply(1:length(l), function(i){
                l2 <- l
                l2[[i]] <- c(l2[[i]], n)
                l2
              })
            }), recursive=FALSE)
          )
        }
      }
    }
  }
  class(m) <- "clusteration"
  m
}

print.clusteration <- function(x, within="-", between=" ", ...){
  f <- function(cell){
    for (i in 1:length(cell))
      cat(i, ": ", paste(sapply(cell[[i]], paste, collapse=within), collapse=between), "\n", sep="")
  }
  if (length(dim(x))==2){
    for (n in 1:nrow(x)){
      for (k in 1:n){
        cat("-- n=", n, " k=", k, " ----------\n", sep="")
        f(x[[n,k]])
      }
    }
  }else{
    n <- attr(x, "n")
    for (k in 1:length(x)){
      cat("-- n=", n, " k=", names(x)[k], " ----------\n", sep="")
      f(x[[k]])
    }
  }
}

