## lotka-volterra code from lecture

run.lv <- function(r,k,a,N0,num.iter) {
  seed <- floor(runif(1)*1e6) ## useful way to track your random seed, for the stochastic version
  N1 <- rep(NA,num.iter)
  N2 <- rep(NA,num.iter)
  N1[1] <- N0[1]
  N2[1] <- N0[2]

  for(t in 2:num.iter) {
    N1[t] <- N1[t-1] * (1+r[1]*(1-(N1[t-1]+a['a12']*N2[t-1])/k[1]))
    N2[t] <- N2[t-1] * (1+r[2]*(1-(N2[t-1]+a['a21']*N1[t-1])/k[2]))
  }
  return(list(N1=N1,N2=N2,seed=seed))
}

## a12>a21
out <- run.lv(r=c(0.3,0.2),
              k=c(100,100),
              a=c(a12=0.2,a21=0.1),
              N0=c(10,20),
              num.iter=5)

## a12<a21
out <- run.lv(r=c(0.3,0.2),
              k=c(100,100),
              a=c(a12=0.2,a21=0.1),
              N0=c(10,20),
              num.iter=5)

plot(out$N1, type='l', col='red', ylim=c(0,max(out$N1,out$N2)))
lines(out$N2, col='blue')

## ------------------------------------------------------------
## code from after class:

## using an apply statement to automate a task
f <- function(x) {
  return(runif(1))
}
f(1)
sapply(1:10, f)
## or, all in one step
sapply(1:10, function(x) runif(1))

## setting the seed:
set.seed(1)
runif(1)

## using parallel apply statements
library(parallel)

## making a slow function
slow.function <- function(x) {
  Sys.sleep(1)
}

## comparing performance on one core
time.start <- proc.time()
out <- lapply(1:16, slow.function)
time.end <- proc.time()
time.end-time.start
## to multiple cores
time.start <- proc.time()
out <- mclapply(1:16, slow.function, mc.cores=8, mc.preschedule=FALSE)
time.end <- proc.time()
time.end-time.start