# comparison of trace maximization as in cl_consensus from R package clue
# with truematch

# R 2.4.1 -- create trace-maximization consensus matrix with clue

library(clue)
R <- 10000
n <- 100

method <- "DWH"  # the default and most similar to what truecluster does (but truecluster has not a small finite ensemble)

e <- do.call("cl_ensemble", lapply(1:R, function(i){
  as.cl_partition(sample(c(rep(1, n-1), 2)))
}))
p1 <- cl_consensus(e, method = method)
p1 <- p1$.Data

e <- do.call("cl_ensemble", lapply(1:R, function(i){
  as.cl_partition(sample(rep(1:2, n/2)))
}))
p2 <- cl_consensus(e, method = method)
p2 <- p2$.Data

e <- do.call("cl_ensemble", lapply(1:R, function(i){
  as.cl_partition(sample(c(rep(1, n/2), rep(2, n/2-1), 3)))
}))
p3 <- cl_consensus(e, method = method)
p3 <- p3$.Data

e <- do.call("cl_ensemble", lapply(1:R, function(i){
  as.cl_partition(c(rep(1, n/2), rep(2, n/2)))
}))
p4 <- cl_consensus(e, method = method)
p4 <- p4$.Data

e <- do.call("cl_ensemble", lapply(1:R, function(i){
  as.cl_partition(c(rep(1, n/2), sample(c(rep(2, n/2-1), 3))))
}))
p5 <- cl_consensus(e, method = method)
p5 <- p5$.Data

save(p1, p2, p3, p4, p5, file="c:/tmp/p.RData")



# R 1.6.2 -- create consensus matrix using truematch
load(file="c:/tmp/p.RData")

library(truecluster)
R <- 10000
n <- 100

#method <- "tracemax"
method <- "truematch"


# NOT using crossprod of soft clustering as in clue
if (FALSE){
  P1 <- matrix(0,n,2)
  P1[cbind(1:n,sample(c(rep(1, n-1), 2)))] <- 1
  for (r in 2:R){
    P <- matrix(0,n,2)
    P[cbind(1:n, sample(c(rep(1, n-1), 2)))] <- 1
    matchind <- matchindex(crossprod(P1,P), method=method)
    P1 <- P1 + P[,matchind$col[matchind$row]]
  }

  P2 <- matrix(0,n,2)
  P2[cbind(1:n,sample(rep(1:2, n/2)))] <- 1
  for (r in 2:R){
    P <- matrix(0,n,2)
    P[cbind(1:n, sample(rep(1:2, n/2)))] <- 1
    matchind <- matchindex(crossprod(P2,P), method=method)
    P2 <- P2 + P[,matchind$col[matchind$row]]
  }

  P3 <- matrix(0,n,3)
  P3[cbind(1:n,sample(c(rep(1, n/2), rep(2, n/2-1), 3)))] <- 1
  for (r in 2:R){
    P <- matrix(0,n,3)
    P[cbind(1:n, sample(c(rep(1, n/2), rep(2, n/2-1), 3)))] <- 1
    matchind <- matchindex(crossprod(P3,P), method=method)
    P3 <- P3 + P[,matchind$col[matchind$row]]
  }

  P4 <- matrix(0,n,3)
  P4[cbind(1:n,c(rep(1, n/2), rep(2, n/2)))] <- 1
  for (r in 2:R){
    P <- matrix(0,n,3)
    P[cbind(1:n, c(rep(1, n/2), rep(2, n/2)))] <- 1
    matchind <- matchindex(crossprod(P4,P), method=method)
    P4 <- P4 + P[,matchind$col[matchind$row]]
  }

  P5 <- matrix(0,n,3)
  P5[cbind(1:n,c(rep(1, n/2), sample(c(rep(2, n/2-1), 3))))] <- 1
  for (r in 2:R){
    P <- matrix(0,n,3)
    P[cbind(1:n, c(rep(1, n/2), sample(c(rep(2, n/2-1), 3))))] <- 1
    matchind <- matchindex(crossprod(P5,P), method=method)
    P5 <- P5 + P[,matchind$col[matchind$row]]
  }
}


# using crosstabulation of hard clustering as in truecluster (gives even more balanced results)
  P1 <- matrix(0,n,2)
  P1[cbind(1:n,sample(c(rep(1, n-1), 2)))] <- 1
  for (r in 2:R){
    c1 <- sample(c(rep(1, n-1), 2))
    i <- cbind(1:n, align(c1, mmcc2cluster(P1)))
    P1[i] <- P1[i] + 1
  }

  P2 <- matrix(0,n,2)
  P2[cbind(1:n,sample(rep(1:2, n/2)))] <- 1
  for (r in 2:R){
    c1 <- sample(rep(1:2, n/2))
    i <- cbind(1:n, align(c1, mmcc2cluster(P2)))
    P2[i] <- P2[i] + 1
  }

  P3 <- matrix(0,n,3)
  P3[cbind(1:n,sample(c(rep(1, n/2), rep(2, n/2-1), 3)))] <- 1
  for (r in 2:R){
    c1 <- sample(c(rep(1, n/2), rep(2, n/2-1), 3))
    i <- cbind(1:n, align(c1, mmcc2cluster(P3)))
    P3[i] <- P3[i] + 1
  }

  P4 <- matrix(0,n,3)
  P4[cbind(1:n,c(rep(1, n/2), rep(2, n/2)))] <- 1
  for (r in 2:R){
    c1 <- c(rep(1, n/2), rep(2, n/2))
    i <- cbind(1:n, align(c1, mmcc2cluster(P4)))
    P4[i] <- P4[i] + 1
  }

  P5 <- matrix(0,n,3)
  P5[cbind(1:n,c(rep(1, n/2), sample(c(rep(2, n/2-1), 3))))] <- 1
  for (r in 2:R){
    c1 <- c(rep(1, n/2), sample(c(rep(2, n/2-1), 3)))
    i <- cbind(1:n, align(c1, mmcc2cluster(P5)))
    P5[i] <- P5[i] + 1
  }

save(p1, p2, p3, p4, p5, P1, P2, P3, P4, P5, file="c:/tmp/p.RData")
load(file="c:/tmp/p.RData")

round(rbind(
  "tracemax random 99:1"=mmcp2ALL(p1)
, "tracemax2 random 50:50"=mmcp2ALL(p2)
, "tracemax3 random 50:49:1"=mmcp2ALL(p3)
, "tracemax4 justified 50:50"=mmcp2ALL(p4)
, "tracemax5 justified 50 random 49:1"=mmcp2ALL(p5)
, "truematch random 99:1"=mmcp2ALL(mmcc2mmcp(P1))
, "truematch2 random 50:50"=mmcp2ALL(mmcc2mmcp(P2))
, "truematch3 random 50:49:1"=mmcp2ALL(mmcc2mmcp(P3))
, "truematch4 justified 50:50"=mmcp2ALL(mmcc2mmcp(P4))
, "truematch5 justified 50 random 49:1"=mmcp2ALL(mmcc2mmcp(P5))
)[,c("H","RMC","I","CIC")], 3)


lapply(1:2, function(i)table(round(p1[,i], 2)))
lapply(1:2, function(i)table(round(mmcc2mmcp(P1)[,i], 2)))
lapply(1:2, function(i)table(round(p2[,i], 1)))
lapply(1:2, function(i)table(round(mmcc2mmcp(P2)[,i], 1)))
lapply(1:3, function(i)table(round(p3[,i], 1)))
lapply(1:3, function(i)table(round(mmcc2mmcp(P3)[,i], 1)))
lapply(1:3, function(i)table(round(p5[1:50,i], 1)))
lapply(1:3, function(i)table(round(p5[51:100,i], 1)))
lapply(1:3, function(i)table(round(mmcc2mmcp(P5)[1:50,i], 1)))
lapply(1:3, function(i)table(round(mmcc2mmcp(P5)[51:100,i], 1)))
