landscape2

library(tree)

## Loading required package: Rcpp

p <- new(Parameters)
p$add_strategy(new(Strategy, list(lma=0.0648406)))
p$add_strategy(new(Strategy, list(lma=0.1977910)))
p$seed_rain <- c(1.1, 1.1)               # Starting rain.
p$seed_rain <- c(348.738300488797, 3.07449195531246) # or start close
p$set_parameters(list(patch_area=1.0))   # See issue #13
p$set_control_parameters(fast.control()) # A bit faster

t.max <- p$disturbance$cdf(tree:::reference.pr.survival.eps)
times0 <- cohort.introduction.times(t.max)
schedule0 <- schedule.from.times(times0, 2L)

# Work out what the equilibrium seed rain is.  This takes a while!  If
# using the second set of seed_rain values this should converge
# quickly once the cohort schedule is worked out.
res <- equilibrium.seed.rain(p, schedule0, 10,
                             build.args=list(verbose=TRUE),
                             progress=TRUE, verbose=TRUE)

## 1: Splitting {75,15} times (141,141)
## 2: Splitting {87,24} times (216,156)
## 3: Splitting {23,16} times (303,180)
## 4: Splitting {10,13} times (326,196)
## 5: Splitting {4,1} times (336,209)
## 6: Splitting {1,0} times (340,210)
## 7: Splitting {1,0} times (341,210)
## 8: Splitting {1,0} times (342,210)
## 9: Splitting {1,0} times (343,210)
## 10: Splitting {1,0} times (344,210)
## *** 1: {348.7,3.074} -> {348.7,3.068} (delta = {-0.0005567,-0.006649})
## *** 2: {348.7,3.068} -> {348.7,3.068} (delta = {-6.16e-05,9.408e-05})
## *** 3: {348.7,3.068} -> {348.7,3.068} (delta = {2.069e-06,-1.329e-06})
## *** 4: {348.7,3.068} -> {348.7,3.068} (delta = {8.311e-08,2.535e-08})
## *** 5: {348.7,3.068} -> {348.7,3.068} (delta = {-2.794e-07,-5.249e-09})
## *** 6: {348.7,3.068} -> {348.7,3.068} (delta = {2.523e-07,1.827e-09})
## *** 7: {348.7,3.068} -> {348.7,3.068} (delta = {8.214e-08,4.248e-09})
## *** 8: {348.7,3.068} -> {348.7,3.068} (delta = {-2.081e-07,-5.33e-09})
## *** 9: {348.7,3.068} -> {348.7,3.068} (delta = {-6.508e-09,1.176e-09})
## *** 10: {348.7,3.068} -> {348.7,3.068} (delta = {-1.5e-07,-1.903e-09})

# Set the seed rain back into the parameters
p$seed_rain <- res$seed_rain[,"out"]

# And rerun again to get the ODE times.
schedule <- res$schedule
ebt <- run.ebt(p, schedule$copy())
schedule$ode_times <- ebt$ode_times

# Resident lma values:
lma.res <- sapply(seq_len(p$size), function(i) p[[i]]$parameters$lma)

OK, good to reasonable accuracy; the seed rain values should be close enough to 1.

cmp <- landscape(lma.res, p, schedule)
cmp - 1 # This should be about zero, plus or minus 1e-5 or so.

## [1]  3.597e-10 -1.726e-05

# Mutant LMA values, in increasing numbers, to test how the time
# requirements scale with the number of strategies.  Should be
# sublinear.
lma.2  <- seq_log(0.03, 0.8, 2)
lma.4  <- seq_log(0.03, 0.8, 4)
lma.8  <- seq_log(0.03, 0.8, 8)
lma.16 <- seq_log(0.03, 0.8, 16)
lma.32 <- seq_log(0.03, 0.8, 32)
lma.64 <- seq_log(0.03, 0.8, 64)

(t.2  <- system.time(w.2  <- landscape(lma.2,  p, schedule)))

##    user  system elapsed 
##  15.779   0.009  15.794

(t.4  <- system.time(w.4  <- landscape(lma.4,  p, schedule)))

##    user  system elapsed 
##  23.952   0.018  23.976

(t.8  <- system.time(w.8  <- landscape(lma.8,  p, schedule)))

##    user  system elapsed 
##  40.582   0.025  40.617

(t.16 <- system.time(w.16 <- landscape(lma.16, p, schedule)))

##    user  system elapsed 
##  76.936   0.095  77.144

(t.32 <- system.time(w.32 <- landscape(lma.32, p, schedule)))

##    user  system elapsed 
## 139.372   0.138 139.569

(t.64 <- system.time(w.64 <- landscape(lma.64, p, schedule)))

##    user  system elapsed 
## 274.087   0.368 274.666

tt <- c(t.2[["elapsed"]],
        t.4[["elapsed"]],
        t.8[["elapsed"]],
        t.16[["elapsed"]],
        t.32[["elapsed"]],
        t.64[["elapsed"]])
n <- 2^seq_along(tt)

plot(n, tt / n, log="x")

plot(w.64 ~ lma.64, type="l", log="x")
abline(v=lma.res, lty=3, col="grey")

# Fitness of the mutants in an empty landscape:
w.empty <- landscape.empty(lma.16, p, schedule)

plot(w.empty ~ lma.16, log="x")