L <- function(priors,logic,cost,temp,n) {

  Ln <- normRow(logic * (1/sum(logic))) # L0

  for (i in 1:n) {
    Sn <- normCol(exp(temp * (log(Ln) - cost))) # Def of Sn
    Ln <- normRow(t(t(Sn) * priors))            # Def of Ln
  }

  list(Ln, Sn) # Returns the (Listener, Speaker) pair

}

st      <- 0.1
mi      <- 0.1
ma      <- 2
lambdas <- seq(mi,ma,st)
costs   <- seq(mi,ma,st)
size    <- length(costs)

Llist <- matrix(rep(NA, size*size), nrow=size)

Llist <- foreach(i=1:size, .combine='rbind') %dopar% {
  foreach(j=1:size, .combine='c') %dopar% {
    list[Ln, Sn] <- L(priors,logic,costs[i]*cost,lambdas[j],iter)
    Ln[1,1]
  }
}

df = data.frame(lambda=rep(lambdas,each=size),
                cost=rep(costs,size),
                p=as.vector(Llist))

ggplot(df,aes(x = lambda, y = cost)) +
  scale_fill_gradientn(colours=cbPalette) +
  geom_raster(aes(fill = p), interpolate = TRUE) +
  coord_equal()

lambdas <- seq(0.05,7.5,0.02)
prior   <- seq(0.05,0.95,0.02)

start_time <- Sys.time()
Llist <- foreach(i=1:length(lambdas), .combine='rbind') %dopar% {
  foreach(j=1:length(prior), .combine='c') %dopar% {
    remain <- (1 - prior[j])/2
    priors <- c(remain, remain, prior[j])
    list[Ln, Sn] <- L(priors,logic,cost,lambdas[i],iter)
    Ln[2,2]
  }
}
end_time <- Sys.time()

df = data.frame(lambda=rep(lambdas,length(prior)),
                prior=rep(prior,each=length(lambdas)),
                p=as.vector(Llist))

ggplot(df,aes(x = lambda, y = prior)) +
  scale_fill_gradientn(colours=cbPalette) +
  geom_raster(aes(fill = p), interpolate = TRUE)

l <- 15
worldNames <- c(seq(1,l))
lex <- matrix(1, l, l)
lex[lower.tri(lex)] <- 0
lexNames   <- worldNames
cost       <- c(log(seq(1,l)))/10
logic      <- t(data.frame(lex))
priors     <- rep((1/length(worldNames)), length(worldNames))
temp       <- 1
iter       <- 10

names(cost)     <- lexNames
colnames(logic) <- worldNames

list[Ln, Sn] <- L(priors,logic,cost,temp,iter)

cost2        <- c("cost", signif(cost, digits=2))
names(cost2) <- c("Lex", lexNames)

L <- function(priors,logic,cost,temp,n) {

  Ln <- normRow(logic * (1/sum(logic))) # L0

  for (i in 1:n) {
    Sn <- normCol(exp(temp * (log(Ln) - cost))) # Def of Sn
    Ln <- normRow(t(t(Sn) * priors))            # Def of Ln
    # This should be replaced with an update that updates all elements of the same equivalence class for a given question.
  }

  list(Ln, Sn) # Returns the (Listener, Speaker) pair

}

# Let's not print too many digits for the human eye
options(digits=2)

normCol = function(m) { # Normalises on columns.
  t(t(m)/colSums(m))
}

normRow = function(m) { # Normalises on rows. This was benchmarked.
  m/rowSums(m)
}

# https://stat.ethz.ch/pipermail/r-help/2004-June/053343.html
# See below for a "why"
list <- structure(NA,class="result")
"[<-.result" <- function(x,...,value) {
   args <- as.list(match.call())
   args <- args[-c(1:2,length(args))]
   length(value) <- length(args)
   for(i in seq(along=args)) {
     a <- args[[i]]
     if(!missing(a)) eval.parent(substitute(a <- v,list(a=a,v=value[[i]])))
   }
   x
}

# Allows returning/matching tuples from functions:
useless   <- function() { list(a = 2, b = 6) }
list[a,]  <- useless() # sets a to 2, ignores "b"
list[,b]  <- useless() # sets b to 6, ignores "a"
list[c,d] <- useless() # sets c to 2 and d to 6