###################################################
require("partitions")
require("gtools")


###################################################
setparts(c(2,2,1))


###################################################
a <- c(1,2,3,2,1)
split(seq_along(a),a)


###################################################
m <- 9
n <- 3
jj <- setparts(restrictedparts(m, n, include.zero=FALSE))
summary(jj)


###################################################
tau <- 1:9
slowest <-  function(x) max(tapply(tau, x, sum))


###################################################
time.taken <- apply(jj, 2, slowest)
minimal.time <- sum(tau)/n
jj[,time.taken == minimal.time]


###################################################
E <- matrix(c(0,0,0,1,1,0,0,0,1,1,1,1,1,0,0), 5,3)
dimnames(E) <-  list(
                     evidence=paste("E",1:5,sep=""),
                     crime=paste("C",1:3,sep="")
                     )
E


###################################################
a <- 
structure(
          c(2, 2, 4, 4, 3, 1, 2, 1, 1, 2, 2, 1, 4, 2, 1, 2, 1, 
            4, 2, 4, 2, 1, 4, 4, 1, 2, 1, 4, 2, 1, 1, 2, 1, 1, 2),
          .Dim = c(5L, 7L),
          .Dimnames = list(
            evidence=paste("E",1:5,sep=""),
            crime=paste("C",1:7,sep="")
            )
)
a


###################################################
genbeta <- function(x){gamma(sum(x))/prod(gamma(x))}

lp <- function(evidence, alpha, partition){
  r <- length(unique(partition))
  evidence <- as.factor(evidence)
  levels(evidence) <- seq_along(alpha)
  out <- 1
  for(k in unique(partition)){  #k is a perp
    thisperp <- partition==k
    evidence.thisperp <- evidence[thisperp]
    no.of.crimes.thisperp <- sum(thisperp)
    out <- out/genbeta(alpha+table(evidence.thisperp))
  }
  return(out*genbeta(alpha)^r)
}

sp <- setparts(7)
l1 <- apply(sp,2, function(x){lp(evidence=a[1,],alpha=rep(1,2),partition=x)})
l2 <- apply(sp,2, function(x){lp(evidence=a[2,],alpha=rep(1,2),partition=x)})
l3 <- apply(sp,2, function(x){lp(evidence=a[3,],alpha=rep(1,4),partition=x)})
l4 <- apply(sp,2, function(x){lp(evidence=a[4,],alpha=rep(1,4),partition=x)})
l5 <- apply(sp,2, function(x){lp(evidence=a[5,],alpha=rep(1,4),partition=x)})

likelihood <- l1*l2*l3*l4*l5
likelihood <- likelihood/max(likelihood)
support <- log(likelihood)


###################################################
par(mfcol = c(1, 2))
plot(likelihood, ylab = "likelihood", xlab = "hypothesis")
abline(h = exp(-2))
plot(support, ylab = "support", xlab = "hypothesis")
abline(h = -2)


###################################################
sp <- setparts(7)
sp[,which.max(support)]

index <- order(support, decreasing = TRUE)
sp <- sp[,index]
support <- support[index]
dimnames(sp) <- list(
                     crime = paste("C",1:7, sep=""),
                     partition =  paste("H",1:ncol(sp),sep="")
                     )

sp[, support > -2]
support[support > -2]