# Written June 22nd, 2006
# Author: Markus Triska triska@metalevel.at
# Public domain code.

# Adjust the matrix m of success probabilities to run your own
# simulations.

m <- matrix(rep(0.5, 32^2), ncol=32)

m[1,] <- 0.9
m[,1] <- 0.1

######################################################################

win <- function(pair) {
  npair <- as.numeric(unlist(strsplit(pair, " ")))
  p1 <- npair[1]
  p2 <- npair[2]
  prob1 <- m[p1,p2]
  w <- sample(c(p1,p2), 1, prob=c(prob1, 1-prob1))
  return(w)
}


winners <- function(group, n) {
  o <- outer(group, group, "paste")
  pairs <- o[upper.tri(o)]
  wins <- unlist(lapply(pairs, win))
  t <- table(wins)
  val <- (as.numeric(names(sort(t, decreasing=T)[1:n])))
  # print(val)
  return(val)
}

worldcup <- function() {

  # stage 1 (group)

  g1 <- winners(1:4, 2)
  g2 <- winners(5:8, 2)

  g3 <- winners(9:12, 2)
  g4 <- winners(13:16, 2)

  g5 <- winners(17:20, 2)
  g6 <- winners(21:24, 2)

  g7 <- winners(25:28, 2)
  g8 <- winners(29:32, 2)

  # stage 2 (knockout)

  k1 <- winners(c(g1[1], g2[2]), 1)
  k2 <- winners(c(g1[2], g2[1]), 1)

  k3 <- winners(c(g3[1], g4[2]), 1)
  k4 <- winners(c(g3[2], g4[1]), 1)

  k5 <- winners(c(g5[1], g6[2]), 1)
  k6 <- winners(c(g5[2], g6[1]), 1)

  k7 <- winners(c(g7[1], g8[2]), 1)
  k8 <- winners(c(g7[2], g8[1]), 1)

  # quarter finals

  q1 <- winners(c(k1, k2), 1)
  q2 <- winners(c(k3, k4), 1)

  q3 <- winners(c(k5, k6), 1)
  q4 <- winners(c(k7, k8), 1)

  # semi finals

  s1 <- winners(c(q1, q2), 1)
  s2 <- winners(c(q3, q4), 1)

  # final

  f <- winners(c(s1, s2), 1)

  return(f)
}

simulations <- function(n) {
  s <- c()
  for (i in 1:n) {
    print(paste("simulating ", i, "/", n))
    s <- c(worldcup(), s)
  }
  return(s)
}