Code I'm not proud of...

I was asked for the code I used to make my little pics. Well, it's not pretty, it's just a few functions in R. For example:

weasel.matrix.solo <- function(target.length=5, alphabet.length=27,mut.prob=.04,mut.effect=FALSE){
if(mut.effect) mut.prob <- alphabet.length*mut.prob/(alphabet.length-1)
M <- matrix(integer((target.length+1)*(target.length+1)),nrow=target.length+1)
for(k in 1:(target.length)){
verschlechtert <- integer(2*target.length+2)
verbessert <- verschlechtert
verschlechtert[1:k + target.length +1] <- dbinom(0:(k-1),k-1,mut.prob*(alphabet.length-1)/alphabet.length)
verbessert[1:(target.length+2-k) + target.length + 1] <- dbinom(0:(target.length+1-k),(target.length+1-k),mut.prob/alphabet.length)
for(kk in 1:(target.length+1)){
M[k,kk] <-sum( verbessert[(0:target.length) + target.length+3-k]*verschlechtert[0:target.length +target.length+3-kk])
}
}
M[target.length+1,target.length+1] <- 1
return(M)
}


This creates a transition matrix for Dawkins's algorithm with a single child. target.length is the length of the target phrase, alphabet.length is self-explaining, mut.prob is the probability with which a letter is mutated: if mut.effect==FALSE, such a mutation has no effect with probability 1/N, i.e., it's a change to the same letter :-)


weasel.matrix.cum <- function(M, pop.size = 10){
target.length <- dim(M)[1]-1
M1 <- M
for(k in 1:(target.length)){
M1[k,] <- cumsum(M[k,])
}
M1 <- M1^pop.size
for(k in 1:(target.length))
M1[k,] <- M1[k,] - c(0,M1[k,-target.length-1])
return(M1)
}

creates a transition matrix for a whole population of size pop.size, given a transition matrix M for a single-child process.

The last thing which is needed is this little snippet:
erw.exakt <- function(M,start.vektor){

M.size <- dim(M)[1]
M <- M[-M.size, -M.size]
M <- solve(diag(M.size-1)-M)
M <- apply(M,1,sum)
return(sum(start.vektor[1:(M.size-1)] * M))
}
Here, the expected value of generations for a transition matrix M is calculated, using a initial distribution for start.vektor (mostly Bin(target.length,mut.prob) distributed - or just the first unit vector).

If I ever prettify this code, I'll take Tom English on his generous offer to put it somewhere on the net.