# R-1.9.1
# Truematch simulations
library(truecluster)
library(e1071)

aapply <- function(dn, FUN, atomic=FALSE, verbose=FALSE, ...)
  # applies function FUN= to each combination of vector elements given in the list of vectors dn=
  # the result is an (list) array of dimension sapply(dn, length)
  # if FUN returns a scalar and atomic=TRUE the result is an atomic array
{
  d <- sapply(dn, length)
  n <- length(d)
  if (atomic){
    a <- array(NA, dim=d, dimnames=dn)
  }else{
    a <- vector(rev(cumprod(d))[1], mode="list")
    dim(a) <- d
    dimnames(a) <- dn
  }
  i <- rep(1, n)
  names(i) <- names(dn)
  repeat{
    arg <- lapply(as.list(1:n), function(j)dn[[j]][[i[[j]]]])
    names(arg) <- names(dn)
    if (verbose)
      cat(paste(names(arg), arg, sep="=", collapse=","), " ... ")
    tim <- system.time(
      x <- do.call("FUN", c(arg, list(...)))
    )[3]
    if (verbose)
      cat(tim, "\n")
    eval(parse(text=paste("a[[", paste(i, collapse=","), "]] <- x", sep="")))
    for (j in 1:n){
      if (i[[j]]==d[[j]]){
        i[[j]] <- 1
      }else{
        i[[j]] <- i[[j]] + 1
        break
      }
    }
    if (all(i==1))
      break
  }
  a
}

matchsim <- function(p, kappa=0.5, R=100, N=100, fixit=TRUE)
{
  p <- c(1-p, p)
  p.0 <- outer(p, p)
  p.1 <- diag(2) * p
  p.cell <- p.1*kappa + p.0*(1-kappa)
  p.cond <- diag(p.cell)/p
  x0 <- rep(factor(1:2), round(p*N))  # this is the true pattern
  f0 <- table(x0)
  cc <- matrix(0, nrow=N, ncol=2)
  cc2 <- matrix(0, nrow=N, ncol=2)
  for (r in 1:R){
    i <- sample(1:N, N, TRUE)
    x <- factor(ifelse(runif(N) < p.cond[x0[i]], x0[i], 3-as.integer(x0[i])), levels=1:2)
    if (fixit){
      f <- table(x)
      d <- f-f0
      if (max(d)>0){
        j <- which.max(d)
        x[x==j][sample(1:f[j], d[j])] <- (2:1)[j]
      }
    }
    if (r>1){
      tab <- table(est[i], x)
      ind <- matchindex(tab, method="heuristic")
      ind2 <- matchindex(tab, method="tracemax")
      x2 <- match(x, ind2$col[match(1:2, ind2$row)])
      x <- match(x, ind$col[match(1:2, ind$row)])
    }else{
      x2 <- x
    }
    x <- x[!duplicated(i)]
    x2 <- x2[!duplicated(i)]
    i <- i[!duplicated(i)]
    cc[cbind(i, x)]  <- cc[cbind(i, x)]  + 1
    cc2[cbind(i, x2)] <- cc2[cbind(i, x2)] + 1
    est <- factor(rowMax(cc, value=FALSE), levels=1:2)
  }
  names(p) <- 1:2
  dimnames(p.cell) <- list(1:2, 1:2)
  ret <- list(true=x0, cc=cc, cc2=cc2, p=p, kappa=kappa, p.cell=p.cell)
  class(ret) <- "matchsim"
  ret
}

summary.matchsim <- function(x, ...){
  p1 <- mmcc2mmcp(x$cc)
  p2 <- mmcc2mmcp(x$cc2)
  groups <- 1:2
  names(groups) <- paste("mean across true cases", groups)
  e1 <- do.call("rbind", lapply(as.list(groups), function(i){
    ret <- colMeans(p1[x$true==i, , drop=FALSE])
    names(ret) <- groups
    ret
  }))
  e2 <- do.call("rbind", lapply(as.list(groups), function(i){
    ret <- colMeans(p2[x$true==i, , drop=FALSE])
    names(ret) <- groups
    ret
  }))
  i1 <- mmcp2I(p1)
  h1 <- mmcp2H(p1)
  cic1 <- i1-h1
  i2 <- mmcp2I(p2)
  h2 <- mmcp2H(p2)
  cic2 <- i2-h2
  ret1 <- c(i=i1, h=h1, cic=cic1)
  ret2 <- c(i=i2, h=h2, cic=cic2)
  rbind("match"=ret1, "match2"=ret2)
}

print.matchsim <- function(x, ...){
  cat("\np\n")
  print(x$p)
  cat("\nf\n")
  print(table(x$true))
  cat("\nh\n")
  print(h(x$p))
  cat("\nkappa\n")
  print(x$kappa)
  cat("\np.cell\n")
  print(x$p.cell)
  p1 <- mmcc2mmcp(x$cc)
  p2 <- mmcc2mmcp(x$cc2)
  groups <- 1:2
  names(groups) <- paste("mean across true cases", groups)
  e1 <- do.call("rbind", lapply(as.list(groups), function(i){
    ret <- colMeans(p1[x$true==i, , drop=FALSE])
    names(ret) <- groups
    ret
  }))
  e2 <- do.call("rbind", lapply(as.list(groups), function(i){
    ret <- colMeans(p2[x$true==i, , drop=FALSE])
    names(ret) <- groups
    ret
  }))
  i1 <- mmcp2I(p1)
  h1 <- mmcp2H(p1)
  cic1 <- i1-h1
  i2 <- mmcp2I(p2)
  h2 <- mmcp2H(p2)
  cic2 <- i2-h2
  ret1 <- c(i=i1, h=h1, cic=cic1)
  ret2 <- c(i=i2, h=h2, cic=cic2)
  ret <- rbind("match"=ret1, "match2"=ret2)
  cat("\nmatch\n")
  print(e1)
  cat("\nmatch2\n")
  print(e2)
  cat("\n\n")
  print(ret)
  invisible(ret)
}

plotit <- function(obj, stat, cont=TRUE, smooth=TRUE, ...){
  oldpar <- par(mfrow=c(2,1), mar=c(0, 4, 1, 0)+0.1, oma=c(3,0,0,0)+0.1)
  on.exit(par(oldpar))
  z <- apply(obj, 1:2, function(e)summary(e[[1]])["match2",stat])
  image(p1, kappa, z, col = rainbow(256), xlab="", ylab="", axes=FALSE)
  axis(2)
  mtext("trace maximization", 3, font=2)
  title(ylab="kappa")
  if (cont){
    if (smooth){
      d <- data.frame(p=as.numeric(rownames(z))[row(z)], kappa=as.numeric(colnames(z))[col(z)], value=as.vector(z))
      m <- loess(value~p+kappa, data=d, ...)
      z2 <- z; z2[] <- predict(m)
      contour(p1, kappa, z2, add=TRUE)
    }else{
      contour(p1, kappa, z, add=TRUE, ...)
    }
  }
  z <- apply(obj, 1:2, function(e)summary(e[[1]])["match",stat])
  image(p1, kappa, z, col = rainbow(256), xlab="", ylab="", axes=FALSE)
  axis(2)
  mtext("truematch algorithm", 3, font=2)
  title(ylab="kappa")
  axis(1, outer=TRUE)
  axis(1, outer=TRUE, at=c(0.12, 0.88), labels=c("p", "p"), tick=FALSE, lty=0)
  if (cont){
    if (smooth){
      d <- data.frame(p=as.numeric(rownames(z))[row(z)], kappa=as.numeric(colnames(z))[col(z)], value=as.vector(z))
      m <- loess(value~p+kappa, data=d, ...)
      z2 <- z; z2[] <- predict(m)
      contour(p1, kappa, z2, add=TRUE)
    }else{
      contour(p1, kappa, z, add=TRUE, ...)
    }
  }
}

sumit <- function(obj, stat){
  m<- apply(obj, 1:2, function(e)summary(e[[1]])["match",stat])
  m2<- apply(obj, 1:2, function(e)summary(e[[1]])["match2",stat])
  list(match=m, match2=m2)
}

# Typical treatment of a non-match
NonMatch <- rbind(
  c(1 , 0)
, c(98, 1)
)
TypicalTracemax  <- matchtable(NonMatch, method="tracemax")
TypicalTruematch <- matchtable(NonMatch)
TypicalTracemax
TypicalTruematch

# Tracemax: high diag fraction and rand index, low kappa and crand
round(unlist(classAgreement(TypicalTracemax)), 6)
# Truematch: low diag fraction, kappa and crand, high rand index
round(unlist(classAgreement(TypicalTruematch)), 6)


# using expectations
RandomMatch <- rbind(
  c(99 , 0)
, c(0, 1)
)
TracemaxNonRandomMatch <- rbind(
  c(98 , 1)
, c(1  , 0)
)
TracemaxExpected <- 0.99 * TracemaxNonRandomMatch + 0.01 * RandomMatch
TracemaxExpected

TruematchNonRandomMatch1 <- rbind(
  c(1 , 0)
, c(98, 1)
)
TruematchNonRandomMatch2 <- rbind(
  c(1 ,98)
, c(0 , 1)
)
TruematchExpected <- 0.99/2 * TruematchNonRandomMatch1 + 0.99/2 * TruematchNonRandomMatch2 + 0.01 * RandomMatch
TruematchExpected

# Tracemax: only random-corrected indices are low (but random correction is inappropriate in cluster context)
round(rbind(
  "TracemaxNonRandomMatch  99%" = unlist(classAgreement(TracemaxNonRandomMatch))
, "           RandomMatch   1%" = unlist(classAgreement(RandomMatch))
, "      TracemaxExpected 100%" = 0.99 * unlist(classAgreement(TracemaxNonRandomMatch)) + 0.01 * unlist(classAgreement(RandomMatch))
), 6)

# Truematch: random-corrected indices AND simple diag fraction is close to zero, only rand is close to 1
round(rbind(
  "TruematchNonRandomMatch1  49.5%" = unlist(classAgreement(TruematchNonRandomMatch1))
, "TruematchNonRandomMatch2  49.5%" = unlist(classAgreement(TruematchNonRandomMatch2))
, "             RandomMatch   1.0%" = unlist(classAgreement(RandomMatch))
, "       TruematchExpected 100.0%" = 0.99/2 * unlist(classAgreement(TruematchNonRandomMatch1)) + 0.99/2 * unlist(classAgreement(TruematchNonRandomMatch2)) + 0.01 * unlist(classAgreement(RandomMatch))
), 6)


p1 <- c(seq(0.01, 0.09, 0.01), seq(0.10, 0.90, 0.02), seq(0.91, 0.99, 0.01))
kappa <- seq(0, 1, 0.025)

# this takes very long
a.nonfix <- aapply(
  list(p=p1, kappa=kappa)
, FUN = matchsim
, verbose = TRUE
, R = 1000
, fixit=FALSE
)
save(a.nonfix, file="matchsimul.RData")
a.fix <- aapply(
  list(p=p1, kappa=kappa)
, FUN = matchsim
, verbose = TRUE
, R = 1000
)
# therefore we save the results before proceeding
save(a.fix, a.nonfix, file="truematch.RData")
load("truematch.RData")

# Generate Figures
postscript("Fig1_Truematch_NonFixI.ps", width=6, height=10)
plotit(a.nonfix, "i", span=0.1)
dev.off()
postscript("Fig2_Truematch_NonFixH.ps", width=6, height=10)
plotit(a.nonfix, "h", span=0.1)
dev.off()
postscript("Fig3_Truematch_NonFixCIC.ps", width=6, height=10)
plotit(a.nonfix, "cic", span=0.1)
dev.off()
postscript("Fig4_Truematch_FixI.ps", width=6, height=10)
plotit(a.fix, "i", span=0.1)
dev.off()
postscript("Fig5_Truematch_FixH.ps", width=6, height=10)
plotit(a.fix, "h", span=0.1)
dev.off()
postscript("Fig6_Truematch_FixCIC.ps", width=6, height=10)
plotit(a.fix, "cic", span=0.1)
dev.off()

# some statistics
sumit(a.fix, "i")
sumit(a.fix, "h")
sumit(a.fix, "cic")
sumit(a.nonfix, "i")
sumit(a.nonfix, "h")
sumit(a.nonfix, "cic")
