rbugs.R

#### configurations  12/15/2003
## BUGS stores the executable of bugs
## workingDir is the directory to save all files, default is ??tempdir()??
## bugsWorkingDir is the directory for wine to use windows type directory
## WINE stores the executable of wine
## useWine = TRUE if use wine

## example usage on linux
## schools.sim <- rbugs(data=schools.data, inits, parameters, "schools.bug", n.chains=3, n.iter=1000, workingDir="/var/tmp/jyan/c/tmp", bugsWorkingDir="c:/tmp", useWine=T, wine="/var/scratch/jyan//wine/wine-20031016/wine", debug=T)


rbugs <- function(data, inits, paramSet, model,
                  ## mcmc options
                  n.chains=1, n.iter=2000, n.burnin=floor(n.iter/2),
                  n.thin=max(1, floor(n.chains*(n.iter-n.burnin)/1000)),
                  ## monitoring options
                  dic=FALSE,
                  ## configuration options
                  debug=FALSE,
                  bugs=Sys.getenv("BUGS"),
                  ##"c:/Program Files/WinBUGS14/WinBUGS14.exe",
                  workingDir = NULL, #getwd(),
                  ##"/var/scratch/jyan/c/tmp", # native
                  bugsWorkingDir, # required argument
                  ##"c:/tmp",
                  useWine = FALSE, 
                  wine = Sys.getenv("WINE"),
                  linbugs = TRUE,
                  cleanBugsWorkingDir = FALSE,
                  genFilesOnly = FALSE,
                  verbose = FALSE,
                  seed=314159
                  ## "/var/scratch/jyan/wine/wine-20031016/wine"
                  ){
  ##  start.time <- Sys.time ()
  os.type <- .Platform$OS.type
  if (os.type == "windows") {
    if (!file.exists(bugs))
      stop(paste("BUGS executable", bugs, "does not exists."))
  }
  else if (os.type == "unix") {
    if (useWine) {
      if (!file.exists(wine))
        stop(paste("wine executable", wine, "does not exists."))
      ## how to check the existence of WinBUGS???
    }
    else { ## use linbugs!
      if (length(bugs) == 0) bugs <- system("which linbugs", TRUE)
      if (length(bugs) == 0)
        stop(paste("BUGS executable", bugs, "does not exists."))
    }
  }
  else warning("This function has not been tested on mac-os.")
  
  ## setup workingDir
  bugsWorkingDir <- filePathAsAbsolute(bugsWorkingDir)
  if (is.null(bugsWorkingDir)) {
    bugsWorkingDir <- tempfile("bugsWorkingDir")
    if (!file.exists(bugwWorkingDir)) dir.create(bugsWorkingDir)
    on.exit(if(cleanBugsWorkingDir) unlink(bugsWorkingDir, TRUE))
  }
  if (is.null(workingDir)) {
    if (useWine) workingDir <- driveTr(bugsWorkingDir, .DriveTable)
    else workingDir <- bugsWorkingDir
  }
  else workingDir <- filePathAsAbsolute(workingDir)
  
  ## prepare the model file by 
  ## making a copy of model to the working directory
  if (!file.exists(model)) stop("Model file doesn't exits.")
  model.file <- file.path(workingDir, "model.txt")
  file.copy(model, model.file, overwrite=TRUE)

  ## prepare the data file
  data.file <- file.path(workingDir, "data.txt")
  genDataFile(data, data.file)

  ## prepare the inits files
  inits.file.stem <- file.path(workingDir, "init")
  genInitsFile(n.chains, inits, inits.file.stem)
  inits.files <- paste(inits.file.stem, 1:n.chains, ".txt", sep="")

  ## prepare the script file
  script.file <- paste(workingDir, "script.txt", sep="/")
  genBugsScript(paramSet, n.chains, n.iter, n.burnin, n.thin, dic,
                model.file, data.file, inits.files,
                workingDir, bugsWorkingDir,
                script.file, debug, useWine, linbugs, seed)

  ## change line breaks from "\n" to "\r\n"
  ## otherwise, linbugs would hang!!
  if (linbugs) {
    ## trLbr(script.file)
    ## spend three hours to figure out that script file doesn't need it!
    trLbr(model.file)
    trLbr(data.file)
    for (i in inits.files) trLbr(i)
  }
  
  ## run bugs
  if (genFilesOnly) {
    cat("Files are generated in", workingDir, "\n")
    return(TRUE)
  }
  if (useWine) script.file <- gsub(workingDir, bugsWorkingDir, script.file)
  runBugs(bugs, script.file, n.chains, workingDir, useWine, wine, linbugs, verbose)

  ## collect the output
  out <- getBugsOutput(n.chains, workingDir, linbugs)
  out
}


genDataFile <- function(dataList, dataFile) {
  if (is.numeric(unlist(dataList))) {
    ## cat(dput2bugs(dataList), file = data.file)
    cat(format4Bugs(dataList), file = dataFile, fill = TRUE)
  }
  else {
    data <- lapply(dataList, get, pos = 1)
    names(data) <- dataList
    ## cat(dput2bugs(data), file = data.file, fill = TRUE)
    cat(format4Bugs(data), file = dataFile, fill = TRUE)
  }
}


genInitsFile <- function(n.chains, inits, initsFileStem) {
  for (i in 1:n.chains) {
    file <- paste(initsFileStem, i, ".txt", sep="")
    if (is.function(inits)) cat(format4Bugs(inits()), file=file, fill = TRUE)
    else cat(format4Bugs(inits[[i]]), file=file, fill = TRUE)
  }
}


## genBugsScript <- function(paramSet,
##                           n.chains,
##                           n.iter,
##                           n.burnin,
##                           n.thin,
##                           dic,
##                           model.file,
##                           data.file,
##                           inits.files,
##                           workingDir=NULL, #getwd(),
##                           ## needs to be readable for BUGS
##                           bugsWorkingDir=getwd(), 
##                           script, #output
##                           debug=FALSE, useWine=FALSE) {
##   if (n.chains != length(inits.files)) stop("length(inits.files) should equal n.chains.")
##   ## n.iter <- n.burnin + n.thin * n.keep

##   ## add deviance to the paramSet list
##   paramSet <- c(paramSet, "deviance")

##   ## setup workingDir 
##   if (is.null(workingDir)) {
##     if (useWine) workingDir <- driveTr(bugsWorkingDir, .DriveTable)
##     else workingDir <- bugsWorkingDir
##   }
##   ## necessary if useWine == TRUE:
##   if (useWine) {
##     model.file <- sub(workingDir, bugsWorkingDir, model.file)
##     data.file <- sub(workingDir, bugsWorkingDir, data.file)
##     for (i in 1:length(inits.files))
##       inits.files[i] <- sub(workingDir, bugsWorkingDir, inits.files[i])
##   }
  
##   ##  history <- paste(bugsWorkingDir, "history.txt", sep="/")
##   coda  <- paste(bugsWorkingDir, "coda", sep="/")
##   logodc <- paste(bugsWorkingDir, "log.odc", sep="/")
##   logfile <- paste(bugsWorkingDir, "log.txt", sep="/")
##   initlist <- paste("inits (", 1:n.chains, ", '", inits.files, "')\n", sep="")
##   savelist <- paste("set (", paramSet, ")\n", sep="")
##   ## write out to script.txt
##   cat (
##        "display ('log')\n",
##        "check ('", model.file, "')\n",
##        "data ('", data.file, "')\n",
##        "compile (", n.chains, ")\n",
##        initlist,
##        "gen.inits()\n",
##        "beg (", ceiling(n.burnin / n.thin) + 1, ")\n",
##        "thin.updater (", n.thin, ")\n",
##        savelist,
##        ## some try update before dic.set()
##        "update (", ceiling(n.burnin / n.thin), ")\n",
##        ifelse(dic, "dic.set()\n", "# dic.set()\n"),
##        "update (", ceiling((n.iter - n.burnin) / n.thin), ")\n",
##        "stats (*)\n",
##        "dic.stats()\n",
##        ## "history (*, '", history, "')\n",
##        "coda (*, '", coda, "')\n",
##        "save ('", logodc, "')\n", 
##        "save ('", logfile, "')\n", file=script, sep="", append=FALSE)
##   if (!debug) cat ("quit ()\n", file=script, append=TRUE)
## }



#### run bugs
runBugs <- function(bugs=Sys.getenv("BUGS"),
                    script,
                    n.chains,
                    workingDir,
                    useWine=FALSE,
                    wine = Sys.getenv("WINE"),
                    linbugs=TRUE,
                    verbose = TRUE) {
#  BUGS <- Sys.getenv("BUGS")
#  if (!file.exists(BUGS)) stop(paste(BUGS, "does not exists."))
  if (!linbugs) {
    if (is.na(pmatch("\"", bugs))) bugs <- paste("\"", bugs, "\"", sep="")
    if (is.na(pmatch("\"", script))) script <- paste("\"", script, "\"", sep="")
    command <- paste(bugs, "/par", script)
  }
  else {
    command <- paste(bugs, "< ", script, "> run.out")
  }
  if (useWine) {
    command <- paste(wine, command)

    ## put a "q" to quit from wine debugg
    q.tmp <- tempfile("q")
    on.exit(unlink(q.tmp))
    cat("q\n", file=q.tmp)
    command <- paste(command, "< ", q.tmp)

    ## redirect the erorr/warning message of Wine
    wine.warn <- tempfile("warn")
    on.exit(unlink(wine.warn))
    command <- paste(command, ">", wine.warn, " 2>&1 ")
  }
  
  
  
  ## clean up previous coda files
  fnames <- getCodaFileNames(n.chains, workingDir, linbugs)
  coda.files <- c(fnames$codaIndexFile, fnames$codaFiles)
##   coda.files <- paste ("coda", 1:n.chains, ".txt", sep="")
##   coda.files <- c("codaIndex.txt", coda.files)
##   coda.files <- file.path(workingDir, coda.files)
  for (i in coda.files) {
    ## cat ("Bugs did not run correctly.\n", file=coda.files[i], append=FALSE)
   if (file.exists(i)) file.remove(i)
  }
  log.file <- file.path(workingDir, "log.txt")
  if (file.exists(log.file)) file.remove(log.file)
  
  ## execute it!
  err <- system(command)
  if (err == -1) stop("System call to BUGS failed.")
  ## show log
  if (verbose) file.show(file.path(workingDir, "log.txt"))

  if (!file.exists(coda.files[1])) 
    stop("BUGS stopped before getting to coda.")
}


#### functions to get the output
getCodaFileNames <- function(n.chains, workingDir, linbugs) {
  CODA <- if (linbugs) "codaCODA" else "coda"
  INDEX <- if (linbugs) "index" else "Index"
  CHAIN <- if (linbugs) "chain" else NULL
  coda  <- file.path(workingDir, CODA)
  codaFiles <- paste(coda, CHAIN, 1:n.chains, ".txt", sep="")
  codaIndexFile <- paste(coda, INDEX, ".txt", sep="")
  list(codaFiles=codaFiles, codaIndexFile=codaIndexFile)
}


getBugsOutput <- function(n.chains, workingDir, linbugs=TRUE) {
  fnames <- getCodaFileNames(n.chains, workingDir, linbugs)
  codaFiles <- fnames$codaFiles
  codaIndexFile <- fnames$codaIndexFile
  codaIndex <- read.table(codaIndexFile, header=FALSE, sep="\t", as.is=TRUE)
  n.keep <- codaIndex[1, 3] - codaIndex[1, 2] + 1
  nodes <- codaIndex[, 1]
  n.param <- length(nodes)
  output <- list()
  for (i in 1:n.chains) {
    foo <- read.table(codaFiles[i], header=FALSE)
    iter <- foo[1:n.keep, 1]
    vals <- matrix(foo[,2], n.keep, n.param)
    dat <- as.data.frame(cbind(iter, vals))
    names(dat) <- c("iter", nodes)
    output[[i]] <- dat
  }    
  output
}