R/jags.R
ca991eae
 jags.model <- function(file, data=sys.frame(sys.parent()), inits,
                        n.chains = 1, n.adapt=1000, nchain)
 {
 
     if (missing(file)) {
         stop("Model file name missing")
     }
     if (!file.exists(file)) {
         stop(paste("Model file \"", file, "\" not found", sep=""))
     }
     if (!missing(nchain)) {
         warning("Argument nchain in jags.model is deprecated. Use n.chains.")
         if (missing(n.chains)) {
             n.chains = nchain
         }
     }
     
     p <- .Call("make_console", PACKAGE="CoGAPS") 
     .Call("check_model", p, file, PACKAGE="CoGAPS")
 
     varnames <- .Call("get_variable_names", p, PACKAGE="CoGAPS")
     if (is.environment(data)) {
         ##Get a list of numeric objects from the supplied environment
         data <- mget(varnames, envir=data, mode="numeric",
                      ifnotfound=list(NULL))
         ##Strip null entries
         data <- data[!sapply(data, is.null)]
     }
     else if (is.list(data)) {
         v <- names(data)
         if (is.null(v)) {
             stop("data must be a named list")
         }
         if (any(nchar(v)==0)) {
             stop("unnamed variables in data list")
         }
         if (any(duplicated(v))) {
             stop("Duplicated names in data list: ",
                  paste(v[duplicated(v)], collapse=" "))
         }
         relevant.variables <- v %in% varnames
         data <- data[relevant.variables]
         unused.variables <- setdiff(v, varnames)
         for (i in seq(along=unused.variables)) {
             warning("Unused variable \"", unused.variables[i], "\" in data")
         }
     }
     else {
         stop("data must be a list or environment")
     }
     
     .Call("compile", p, data, as.integer(n.chains), TRUE, PACKAGE="CoGAPS")
 
 ### Setting initial values
 
     if (!missing(inits)) {
 
         checkParameters <- function(inits) {
             if(!is.list(inits))
                 return (FALSE)
 
             inames <- names(inits)
             if (is.null(inames) || any(nchar(inames) == 0))
                 return (FALSE)
 
             if (any(duplicated(inames)))
                 return (FALSE)
             
             if (any(inames==".RNG.name")) {
                 rngname <- inits[[".RNG.name"]]
                 if (!is.character(rngname) || length(rngname) != 1)
                     return (FALSE)
                 inits[[".RNG.name"]] <- NULL
             }
 
             if (!all(sapply(inits, is.numeric)))
                 return (FALSE)
             
             return (TRUE)
         }
         
         setParameters <- function(inits, chain) {
             if (!is.null(inits[[".RNG.name"]])) {
                 .Call("set_rng_name", p, inits[[".RNG.name"]],
                       as.integer(chain), PACKAGE="CoGAPS")
                 inits[[".RNG.name"]] <- NULL
             }
             .Call("set_parameters", p, inits, as.integer(chain),
                   PACKAGE="CoGAPS")
         }
         
         init.values <- vector("list", n.chains)
         
         if (is.function(inits)) {
             if (any(names(formals(inits)) == "chain")) {
                 for (i in 1:n.chains) {
                     init.values[[i]] <- inits(chain=i)
                 }
             }
             else {
                 for (i in 1:n.chains) {
                     init.values[[i]] <- inits()
                 }
             }
         }
         else if (is.list(inits)) {
 
             if (checkParameters(inits)) {
                 ## Replicate initial values for all chains
                 for (i in 1:n.chains) {
                     init.values[[i]] <- inits
                 }
             }
             else {
                 if (length(inits) != n.chains) {
                     stop("Length mismatch in inits")
                 }
                 init.values <- inits
             }
         }
             
         for (i in 1:n.chains) {
             if (!checkParameters(init.values[[i]])) {
                 stop("Invalid parameters for chain ", i)
             }
             setParameters(init.values[[i]], i)
         }
     }
 
     .Call("initialize", p, PACKAGE="CoGAPS")
 
     model.state <- .Call("get_state", p, PACKAGE="CoGAPS")
     model.data <- .Call("get_data", p, PACKAGE="CoGAPS")
     model.code <- readLines(file, warn=FALSE)
     model <- list("ptr" = function() {p},
                   "data" = function() {model.data},
                   "model" = function() {model.code},
                   "state" = function(internal=FALSE)
                   {
                       if(!internal) {
                           for(i in 1:n.chains) {
                               model.state[[i]][[".RNG.state"]] <- NULL
                               model.state[[i]][[".RNG.name"]] <- NULL
                           }
                       }
                       return(model.state)
                   },
                   "nchain" = function()
                   {
                       .Call("get_nchain", p, PACKAGE="CoGAPS")
                   },
                   "iter" = function()
                   {
                       .Call("get_iter", p, PACKAGE="CoGAPS")
                   },
                   "sync" = function() {
                       
                       model.state <<- .Call("get_state", p, PACKAGE="CoGAPS")
                   },
                   "recompile" = function() {
                       ## Clear the console
                       .Call("clear_console", p, PACKAGE="CoGAPS")
                       p <<- .Call("make_console", PACKAGE="CoGAPS")
                       ## Write the model to a temporary file so we can re-read it
                       mf <- tempfile()
                       writeLines(model.code, mf)
                       .Call("check_model", p, mf, PACKAGE="CoGAPS")
                       unlink(mf)
                       ## Re-compile
                       .Call("compile", p, data, n.chains, FALSE, PACKAGE="CoGAPS")
                       ## Re-initialize
                       if (!is.null(model.state)) {
                           if (length(model.state) != n.chains) {
                               stop("Incorrect number of chains in saved state")
                           }
                           for (i in 1:n.chains) {
                               statei <- model.state[[i]]
                               rng <- statei[[".RNG.name"]]
                               if (!is.null(rng)) {
                                   .Call("set_rng_name", p, rng, i, PACKAGE="CoGAPS")
                                   statei[[".RNG.name"]] <- NULL
                               }
                               .Call("set_parameters", p, statei, i, PACKAGE="CoGAPS")
                           }
                           .Call("initialize", p, PACKAGE="CoGAPS")
                           ## Redo adaptation
                           adapting <- .Call("is_adapting", p, PACKAGE="CoGAPS")
                           if(n.adapt > 0 && adapting) {
                               cat("Adapting\n")
                               .Call("update", p, n.adapt, PACKAGE="CoGAPS")
                               if (!.Call("adapt_off", p, PACKAGE="CoGAPS")) {
                                   warning("Adaptation incomplete");
                               }
                           }
                           model.state <<- .Call("get_state", p, PACKAGE="CoGAPS")
                       }
                       invisible(NULL)
                   })
     class(model) <- "jags"
 
     if (n.adapt > 0) {
         adapt(model, n.adapt)
     }
     return(model)
 }
 
 parse.varname <- function(varname) {
 
   ## Try to parse string of form "a" or "a[n,p:q,r]" where "a" is a
   ## variable name and n,p,q,r are integers
 
   v <- try(parse(text=varname, n=1), silent=TRUE)
   if (!is.expression(v) || length(v) != 1)
     return(NULL)
 
   v <- v[[1]]
   if (is.name(v)) {
     ##Full node array requested
     return(list(name=deparse(v)))
   }
   else if (is.call(v) && identical(deparse(v[[1]]), "[") && length(v) > 2) {
     ##Subset requested
     ndim <- length(v) - 2
     lower <- upper <- numeric(ndim)
     if (any(nchar(sapply(v, deparse)) == 0)) {
       ##We have to catch empty indices here or they will cause trouble
       ##below
       return(NULL)
     }
     for (i in 1:ndim) {
       index <- v[[i+2]]
       if (is.numeric(index)) {
         ##Single index
         lower[i] <- upper[i] <- index
       }
       else if (is.call(index) && length(index) == 3 &&
                identical(deparse(index[[1]]), ":") &&
                is.numeric(index[[2]]) && is.numeric(index[[3]]))
         {
           ##Index range
           lower[i] <- index[[2]]
           upper[i] <- index[[3]]
         }
       else return(NULL)
     }
     if (any(upper < lower))
       return (NULL)
     return(list(name = deparse(v[[2]]), lower=lower, upper=upper))
   }
   return(NULL)
 }
 
 parse.varnames <- function(varnames)
 {
   names <- character(length(varnames))
   lower <- upper <- vector("list", length(varnames))
   for (i in seq(along=varnames)) {
     y <- parse.varname(varnames[i])
     if (is.null(y)) {
       stop(paste("Invalid variable subset", varnames[i]))
     }
     names[i] <- y$name
     if (!is.null(y$lower)) {
       lower[[i]] <- y$lower
     }
     if (!is.null(y$upper)) {
       upper[[i]] <- y$upper
     }
   }
   return(list(names=names, lower=lower, upper=upper))
 }
 
 
 jags.samples <-
   function(model, variable.names, n.iter, thin=1, type="trace", ...)
 {
     if (class(model) != "jags")
       stop("Invalid JAGS model")
 
     if (!is.character(variable.names) || length(variable.names) == 0)
       stop("variable.names must be a character vector")
      
     if (!is.numeric(n.iter) || length(n.iter) != 1 || n.iter <= 0)
       stop("n.iter must be a positive integer")
     if (!is.character(type))
       stop("type must be a character vector")
 
     pn <- parse.varnames(variable.names)
     .Call("set_monitors", model$ptr(), pn$names, pn$lower, pn$upper,
           as.integer(thin), type, PACKAGE="CoGAPS")
     update(model, n.iter, ...)
     ans <- .Call("get_monitored_values", model$ptr(), type, PACKAGE="CoGAPS")
     for (i in seq(along=variable.names)) {
       .Call("clear_monitor", model$ptr(), pn$names[i], pn$lower[[i]],
             pn$upper[[i]], type, PACKAGE="CoGAPS")
     }
     return(ans)
 }
 
 nchain <- function(model)
 {
     if (!inherits(model, "jags"))
       stop("Invalid JAGS model object in nchain")
     
     .Call("get_nchain", model$ptr(), PACKAGE="CoGAPS")
 }
 
 load.module <- function(name, path, quiet=FALSE)
 {
     if (name %in% list.modules()) {
         ## This is a stop-gap measure as JAGS 2.1.0 does allow you
         ## to load the same module twice. This should be fixed in
         ## later versions.
         return(invisible()) #Module already loaded
     }
     
     if (missing(path)) {
         path = getOption("jags.moddir")
         if (is.null(path)) {
             stop("option jags.moddir is not set")
         }
     }
     if (!is.character(path) || length(path) != 1)
         stop("invalid path")
     if (!is.character(name) || length(name) != 1)
         stop("invalid name")
 
     file <- file.path(path, paste(name, .Platform$dynlib.ext, sep=""))
     if (!file.exists(file)) {
         stop("File not found: ", file)
     }
     if (!isDLLLoaded(file)) {
         ## We must avoid calling dyn.load twice on the same DLL This
         ## may result in the DLL being unloaded and then reloaded,
         ## which will invalidate pointers to the distributions,
         ## functions and factories in the module.
         dyn.load(file)
     }
     ok <- .Call("load_module", name, PACKAGE="CoGAPS")
     if (!ok) {
         stop("module", name, "not found\n", sep=" ")
     }
     else if (!quiet) {
         cat("module", name, "loaded\n", sep=" ")
     }
     invisible()
 }
 
 list.modules <- function()
 {
     .Call("get_modules", PACKAGE="CoGAPS");
 }
 
 isDLLLoaded <- function(file)
 {
     dll.list <- getLoadedDLLs()
     for (i in seq(along=dll.list)) {
         if (dll.list[[i]]["path"][1] == file)
             return(TRUE)
     }
     return(FALSE)
 }