# Get the lecuyerseed object
set.seed(123321)
nrep <- 1000
library(parallel)
RNGkind("L'Ecuyer-CMRG") # Use L'Ecuyer method
s <- .Random.seed # Extract the current seed number
lecuyerseed <- NULL
for (i in 1:nrep) {
  lecuyerseed <- rbind(lecuyerseed, s) # Save seed number to the list
  s <- nextRNGStream(s) # Build independent seed number
}
save(lecuyerseed, file = "lecuyerseed.Rda")


# Create a function to run an analysis for one condition
runrep <- function(cond, seednum) {
  RNGkind("L'Ecuyer-CMRG")
  assign(".Random.seed", seednum[unlist(cond[4]),], envir = .GlobalEnv)
  a <- unlist(cond[1])
  b <- unlist(cond[2])
  n <- unlist(cond[3])
  x <- rnorm(n, 0, 1)
  pred <- a + b * x
  predprob <- 1/(1+exp(-pred))
  y <- runif(length(predprob)) < predprob
  dat <- data.frame(x, y)
  mylogit <- glm(y ~ x, data = dat, family = "binomial")
  result <- matrix(NA, 3, 2)
  result[1,] <- try(confint.default(mylogit)[2,])
  result[2,] <- try(confint(mylogit)[2,])
  confint.boot <- function(mydata) {
    FUN <- function(data, i) {
      mylogit <- glm(y ~ x, data = data[i,], family = "binomial")
      coef(mylogit)[2]
    }
    library(boot)
    bootout <- boot(data = mydata, statistic = FUN, R = 2000, stype = "i")
    boot.ci(boot.out = bootout, conf = 0.95, type = "bca")[["bca"]][1, 4:5]
  }
  result[3,] <- try(confint.boot(dat))
  cover <- (result[,1] < b) & (b < result[,2])
  return(cover)
}

# Try to run one condition
RNGkind("L'Ecuyer-CMRG")
runrep(c(-2, 1, 200, 50), seednum = lecuyerseed)


# Try to run all conditions
a <- c(-10, -2, 0, 2, 10)
b <- c(0, 1, 2, 3)
n <- c(50, 100, 200, 400, 800)
rep <- 50
conds <- expand.grid(a, b, n, rep)
conds <- data.frame(t(conds))
result <- lapply(conds, runrep, seednum = lecuyerseed)
result <- do.call(rbind, result)
save(result, file=paste("result50.rda", sep=""))

# Create the distributeFile function
distributeFile <- function(startNum, endNum) {

file <- '
# The built package is required to get lecuyerseed. 
library(mypackage) 
library(parallel)
a <- c(-10, -2, 0, 2, 10)
b <- c(0, 1, 2, 3)
n <- c(50, 100, 200, 400, 800)
rep <- subRep
conds <- expand.grid(a, b, n, rep)
conds <- data.frame(t(conds))
result <- mclapply(conds, runrep, seednum = lecuyerseed, mc.cores = 11)
result <- do.call(rbind, result)
save(result, file=paste("result", subRep, ".rda", sep=""))
'

for(i in startNum:endNum) {
  file2 <- gsub("subRep", i, file)
  write(file2, file=paste("sim", i, ".R", sep=""))
}
}

# Start creating a package
package.skeleton(list = c("distributeFile", "runrep", "lecuyerseed"), 
	name = "mypackage")
	
# Install the package
install.packages("mypackage_1.0.tar.gz", repos = NULL, type = "source")