par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0))

model <- function(t, state, parms){
  state <- ifelse(state < 0, 0, state)
  R <- state[1:nr]
  if (nn == 0) return(list(D*(S-R)))
  N <- state[(nr+1):(nr+nn)]
  mu <- r*sapply(seq(nn),function(i){min(R/(Ks[[i]]+R))})
  co <- sapply(seq(nn),function(i){Cs[[i]]*mu[i]*N[i]})
  dR <- D*(S-R) - rowSums(co)
  dN <- 1e-3 + (mu - D)*N
  return(list(c(dR,dN)))    
} 

Ninvaders <- 10
nr <- 3; nn <- 0
p <- NULL
D <- 0.25; Kmean <- 0.5; Smean <- 10; rmean <- 1; CK <- 20
S <- rnorm(nr,Smean,Smean/10)
R <- S; names(R) <- paste("R",seq(1,nr),sep="")
s <- R
run(main=0)

r <- c(); Ks <- list(); Cs <- list()
nreject <- 0; Rrejected <- c(); Krejected <- list(); 

for (i in seq(Ninvaders)) {
  invader <- FALSE
  while (!invader) {
    ri <- rnorm(1,rmean,rmean/10)
    Ki <- runif(nr,0,1)
    Ki <- Kmean*Ki/mean(Ki)
    R0 <- ri*min(s[1:nr]/(Ki+s[1:nr]))/D
    if (R0 > 1) invader <- TRUE
    else {
      nreject <- nreject + 1
      Rrejected <- c(Rrejected,r)
      Krejected[[nreject]] <- Ki
    }
  }
  r <- c(r, ri)
  Ks[[i]] <- Ki
  Cs[[i]] <- Ki*rnorm(nr,1,0.01)/CK
  #Cs[[i]] <- runif(nr,0.04,0.06)
  s <- c(s,0.1)
  names(s) <- c(names(s[1:(nr+nn)]),paste("N",i,sep=""))
  nn <- nn + 1
  s <- run(1e4,main=i)
  print(c(Fitness=R0,s[s>0.1]))
}

present <- which(s[(nr+1):(nr+nn)] > 0.1)
print(c(Ninvasions=i,Npresent=length(present),Nrejected=nreject))
print(c(meanR=mean(r),meanRpresent=mean(r[present]),meanRrej=mean(Rrejected)))
Kpresent <- Ks[present]
print(c(sdK=sd(unlist(Ks)),sdKpresent=sd(unlist(Kpresent)),sdKrejected=sd(unlist(Krejected))))