# truecluster/R/matchindex.R
# Truecluster matching algorithms and heuristics
# (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
###########################################################################

if (FALSE){

  matchindex.R.recall <- function(observed)
    # matches two cluster solutions using a recursive algorithm based on chisquare evaluation
    # observed= input matrix representing cross table of cluster memberships
    # is designed to work with square tables but will also help to match in non-square
    # returns cbind(rowindex, colindex) where the row sequence gives the matching sequence
  {
    rown <- dim(observed)[1]
    coln <- dim(observed)[2]
    if (min(rown, coln)>1) {
      rsum <- observed%*%rep(1, coln)
      csum <- rep(1,rown) %*% observed
      ## used ChiSquare logic: signed standardized residuals
      expected <- (rsum %*% csum) / sum(rsum)
      ominuse <- observed - expected
      signed.standardized.residual <- (ominuse)^2/expected * sign(ominuse)
      i <- order(signed.standardized.residual, observed, runif(rown*coln))[length(observed)]
      rowi <- row(observed)[i]
      coli <- col(observed)[i]
      rest <- Recall(observed[-rowi,-coli, drop=FALSE])
      return(list(row=c(rowi, (1:rown)[-rowi][rest$row]), col=c(coli, (1:coln)[-coli][rest$col])))
    } else {
      return(list(row=1:rown, col=1:coln))
    }
  }

  matchindex.R.while <- function(observed)
  {
    d <- dim(observed)
    rn1 <- rn <- d[1]
    cn1 <- cn <- d[2]
    rpos <- sample(seq(length=rn))
    cpos <- sample(seq(length=cn))
    while(rn>1 && cn>1){
      rsum <- observed[rpos[1:rn], cpos[1:cn]] %*% rep(1, cn)
      csum <- (rep(1,rn) %*% observed[rpos[1:rn], cpos[1:cn]]) / sum(rsum)
      expected <- (rsum %*% csum)
      res <- observed[rpos[1:rn], cpos[1:cn]] - expected
      ssres <- (res)^2 / expected * sign(res)
      #if (rn1==rn){
      #  ssres1 <- ssres[match(1:rn, rpos), match(1:cn, cpos)]
      #}
      #i <- order(ssres1[rpos[1:rn], cpos[1:cn]], observed[rpos[1:rn], cpos[1:cn]], runif(rn*cn))[rn*cn]
      #plot(ssres, ssres1[rpos[1:rn], cpos[1:cn]])
      #abline(0,1)
      i <- order(ssres, observed[rpos[1:rn], cpos[1:cn]], runif(rn*cn))[rn*cn]
      ri <- row(observed[rpos[1:rn], cpos[1:cn]])[i]
      ci <- col(observed[rpos[1:rn], cpos[1:cn]])[i]
      tmp <- rpos[rn]
      rpos[rn] <- rpos[ri]
      rpos[ri] <- tmp
      tmp <- cpos[cn]
      cpos[cn] <- cpos[ci]
      cpos[ci] <- tmp
      rn <- rn - 1
      cn <- cn - 1
    }
    list(row=rev(rpos), col=rev(cpos))
  }

}


#! \name{matchindex}
#! \alias{matchindex}
#! \alias{matchtable}
#! \alias{matchcalcs}
#! \title{ Matching for cross-tables }
#! \description{
#!   Permutes rows and columns of a cross-table to match the marginal classes.
#! }
#! \usage{
#! matchindex(observed, method = c("heuristic", "truematch", "tracemax")[1])
#! matchtable(observed, method = c("heuristic", "truematch", "tracemax")[1])
#! matchcalcs(k)
#! }
#! \arguments{
#!   \item{observed}{ an integer matrix representing a cross-table, non-square allowed }
#!   \item{method}{ one of c("heuristic","truematch","tracemax"), see details }
#!   \item{k}{ number of categories to match }
#! }
#! \details{
#!   Method \emph{heuristic} -- the default -- matches using the truematch heuristic, i.e. recursively calculating resisuals and removing the row/col with the biggest residuals. The number of residuals calculated during matching is given by function \code{matchcalcs}  = 1/6*(n*(n+1)*(2*n+1)) - 1 (Otto, Forster, Analysis 1, Diffential und Integralrechnung einer Veränderlichen, Viehweg, Braunschweig (1992).). Thus the truematch heuristic has time-complexity O(n^3), and memory-complexity O(n^2), if replacing recursion by while loops as done here.
#!   Method \emph{truematch} -- the truematch algorithm -- calculates residuals just once and then applies trace maximization (Munkres Hungarian Method).
#!   Method \emph{tracemax} applies trace maximization directly without calculating residuals.
#!   All methods can cope with (non-square) rectangular matrices. All methods break ties by random. All methods use C implementations. All methods have polynomial time complexity.
#! }
#! \value{
#!   Function \code{matchindex} returns a list with components
#!   \item{row}{permutation index for the rows}
#!   \item{col}{permutation index for the columns}
#!   Function \code{matchtable} returns reordered table
#!   \item{row}{permutation index for the rows}
#!   \item{col}{permutation index for the columns}
#!   Function \code{matchcalcs} returns the number of residuals to calculate during matching
#! }
#! \references{ xx }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{munkres}}, for another implementation  of the Hungarian method see \code{\link[clue]{solve_LSAP}}, for greedy heuristics see \code{\link[e1071]{matchClasses}} }
#!
#! \examples{
#! m <- rbind(c(1, 98), c(0, 1))
#! i <- matchindex(m)
#! i
#! m[i$row, i$col]
#!i <- matchindex(m, method="truematch")
#! i
#! m[i$row, i$col]
#! i <- matchindex(m, method="tracemax")
#! i
#! m[i$row, i$col]
#!   \dontrun{
#!      library(e1071)
#!      library(clue)
#!      #library(lpSolve)
#!
#!      cat("Check quality with respect to trace maximization\n")
#!      x <- integer(6)
#!      names(x) <- c("e","s","l","t","r","m")
#!      for (i in 1:400){
#!         k <- 5
#!         m <- table(sample(1:k, 10000, TRUE), sample(1:k, 10000, TRUE))
#!         i.e <- matchClasses(m, method="exact")
#!         m.e <- m[,i.e]
#!         i.s <- solve_LSAP(m, maximum=TRUE)
#!         m.s <- m[,i.s]
#!         i.l <- apply(lp.assign(-m)$solution, 1, which.max)
#!         m.l <- m[,i.l]
#!         i.t <- matchindex(m, method="tracemax")
#!         m.t <- m[i.t$row, i.t$col]
#!         i.r <- matchindex(m, method="truematch")
#!         m.r <- m[i.r$row, i.r$col]
#!         i.m <- matchindex(m, method="truecluster")
#!         m.m <- m[i.m$row, i.m$col]
#!         x["e"] <- x["e"] + sum(diag(m.e))
#!         x["s"] <- x["s"] + sum(diag(m.s))
#!         x["l"] <- x["l"] + sum(diag(m.l))
#!         x["t"] <- x["t"] + sum(diag(m.t))
#!         x["r"] <- x["r"] + sum(diag(m.r))
#!         x["m"] <- x["m"] + sum(diag(m.m))
#!      }
#!      x / max(x)
#!
#!      cat("Check speed\n")
#!      K <- 2^(2:9)
#!      n <- length(K)
#!      tim <- matrix(NA, nrow=n, ncol=3, dimnames=list(NULL, c("s","m","t")))
#!      for (i in 1:n){
#!         k <- K[i]
#!         m <- table(sample(1:k, 1000000, TRUE), sample(1:k, 1000000, TRUE))
#!         tim[i,"s"] <- system.time({solve_LSAP(m)})[3]
#!         tim[i,"m"] <- system.time({matchindex(m, method="truecluster")})[3]
#!         tim[i,"t"] <- system.time({matchindex(m, method="tracemax")})[3]
#!         matplot(K, tim, pch=c("s","m","t"))
#!      }
#!   }
#! }
#! \keyword{ utilities }

matchindex <- function(observed, method=c("heuristic","truematch","tracemax")[1])
{
  method <- match.arg(method, c("heuristic","truematch","tracemax"))
  if (any(is.na(observed)))
    stop("NAs not allowed")
  d <- dim(observed)
  if (method=="heuristic"){
    ret <- .C("R_tc_matchindex"
    , observed = as.integer(observed)
    , rn   = as.integer(d[1])
    , cn   = as.integer(d[2])
    , rpos = as.integer(sample(seq(length=d[1]))-1)
    , cpos = as.integer(sample(seq(length=d[2]))-1)
    , NAOK = FALSE, DUP = FALSE, PACKAGE = "truecluster"
    )[c("rpos", "cpos")]
    list(row=rev(ret$rpos)+1, col=rev(ret$cpos)+1)
  }else if (method=="truematch"){
    rsum <- observed%*%rep(1, d[2])
    csum <- rep(1,d[1]) %*% observed
    expected <- (rsum %*% csum) / sum(rsum)
    ominuse <- observed - expected
    signed.standardized.residual <- (ominuse)^2/expected * sign(ominuse)
    munkres(-signed.standardized.residual)
  }else{
    observed <- as.double(observed)
    dim(observed) <- d
    munkres(-observed)
  }
}

matchtable <- function(observed, method=c("heuristic","truematch","tracemax")[1])
{
  ind <- matchindex(observed, method=method)
  observed[ind$row, ind$col]
}

matchcalcs <- function(k){
  ((k*(k+1))*(2*k+1))/6-1
}



#! \name{align}
#! \alias{align}
#! \alias{align.default}
#! \alias{compare}
#! \alias{compare.default}
#! \title{ Cluster aligning/comparing }
#! \description{
#!   Aligns or compares objects that are clusterings or contain clusterings
#! }
#! \synopsis{
#!  align(obj, \dots)
#!  compare(obj, \dots)
#! }
#! \usage{
#! align.default(obj, refobj, method = c("heuristic", "truematch", "tracemax")[1], \dots)
#! compare.default(obj, method = c("heuristic", "truematch", "tracemax")[1], \dots)
#! }
#! \arguments{
#!   \item{obj}{ object being or containing a clustering }
#!   \item{refobj}{ an object which defines how to align }
#!   \item{obj2}{ an object to compare with (columns of table) }
#!   \item{method}{ one of c("heuristic","truematch","tracemax"), see \code{\link{matchindex}} }
#!   \item{\dots}{ further arguments depending on class method }
#! }
#! \value{
#!   Function \code{align} returns an aligned object (clustering or object containing a clustering)
#!   Function \code{compare} returns an aligned table
#! }
#! \author{ Jens Oehlschlägel }
#! \seealso{ \code{\link{matchindex}} }
#!
#! \examples{
#!  x <- c(1,1,2,2,3,3)
#!  y <- c(2,2,3,3,1,1)
#!  table(x, y)
#!  align(y, x)
#!  table(x, align(y,x))
#!  compare(x, y)
#! }
#! \keyword{ utilities }
#! \keyword{ cluster }


align <- function(
  obj
, ...
)
UseMethod("align")

compare <- function(
  obj
, ...
)
UseMethod("compare")

align.default <- function(obj, refobj, method=c("heuristic","truematch","tracemax")[1], ...)
  # convenience function: RETURNs cluster aligned with referencecluster
{
      matchtab <- table(refobj, obj, dnn=NULL)
      matchind <- matchindex(matchtab, method=method)
      match(obj, matchind$col[match(1:length(matchind$row), matchind$row)])
}

compare.default <- function(obj, obj2, method=c("heuristic","truematch","tracemax")[1], ...)
  # convenience function: might change
{
      matchtable(table(obj, obj2, dnn=NULL), method=method)
}
