axTexpr <- function(side, at = axTicks(side, axp=axp, usr=usr, log=log), axp = NULL, usr = NULL, log = NULL) { ## Purpose: Do "a 10^k" labeling instead of "a e" ## this auxiliary should return 'at' and 'label' (expression) ## ---------------------------------------------------------------------- ## Arguments: as for axTicks() ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 7 May 2004, 18:01 eT <- floor((abs(at))) mT <- at / 10^eT ss <- lapply(seq(along = at), function(i) substitute(10^E, list(E=eT[i]))) do.call("expression", ss) } N <- 1000 R <- 1000 p <- 0.5; q <- 1 - p K <- 4 * p * q * N * R / (2* N + R) u <- 1 - 1 / (2 * N) - 1 / R t <- 0:5000 J <- K - K *(1 - 2 * p * q / K) ^ t plot(J ~ t, type = "l", xlab = "Generations", ylab = "Number of Junctions") par(new = T) vx <- log(u ^ t - 1 / K) / (log(u) * t) - 1 plot(vx ~ t, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "", lty = 2, col = "blue", ylim = c(0, 1)) axis(4) mtext("Relative Error", side = 4, line = 2.5) R <- 1000 N <- 10^seq(1, 7, by = 0.1) p <- 0.5; q <- 1-p K <- 4*p*q*N*R/(2*N+R) u <- 1 - 1/(2*N) - 1/R MAT <- log(1/K) / log(u) plot(MAT~log10(N), type = "l", ylab = "Maximum Accurate Time", xlab = "Population Size (N)", xaxt = "n") aY <- axTicks(1) aY <- aY[which(floor((aY)) ==(aY))] axis(1, at = aY, label = axTexpr(2, aY), las = 2) par(mfrow=c(1, 1)) N <- 1000 R <- 10^seq(1, 7, by = 0.1) p <- 0.5; q <- 1-p K <- 4*p*q*N*R/(2*N+R) u <- 1 - 1/(2*N) - 1/R MAT <- log(1/K) / log(u) plot(MAT~log10(R), type = "l", ylab = "Maximum Accurate Time", xlab = "Chromosome Size (R)", xaxt = "n") aY <- axTicks(1) aY <- aY[which(floor((aY)) == (aY))] axis(1, at = aY, label = axTexpr(2, aY), las = 2) library(akima) library(fields) r <- 10^seq(1, 10, by = 0.1) n <- 10^seq(1, 10, by = 0.1) p <- 0.5; q <- 1-p; found <- c(); for(N in n) { for(R in r) { u <- 1 - 1/(2*N) - 1/R K <- 4*p*q*N*R/(2*N+R) MAT <- log(1/K) / log(u) if(!is.nan(MAT) && !is.na(MAT) && MAT > 0) { toAdd <- c(N,R,MAT) cat(toAdd,"\n") found <- rbind(found,toAdd); } } } grid <- interp(x = log10(found[, 1]), y = log10(found[, 2]), z = log10(found[, 3]), duplicate = "mean") filled.contour(grid, nlevels = 10, plot.axes = { contour( grid, nlevels = 10, drawlabels = F, axes = FALSE, frame.plot = FALSE, add = TRUE); axis(1, at = aY, label= axTexpr(2, aY), las = 0); axis(2, at = aX, label= axTexpr(2, aX), las = 2) }, col = tim.colors(20), xlab = "Population Size (N)", ylab = "Chromosome Size (R)")