From 2f5a34da83184f014723f6d8cd06a6f824c5b41e Mon Sep 17 00:00:00 2001 From: Jeffrey Pullin Date: Fri, 29 Jan 2021 12:44:12 +1100 Subject: [PATCH] Remove rstan3 directory --- rstan3/.Rbuildignore | 23 -- rstan3/DESCRIPTION | 26 -- rstan3/NAMESPACE | 12 - rstan3/R/AllClass.R | 450 ---------------------------------- rstan3/R/StanFitMCMCMethods.R | 93 ------- rstan3/R/StanTypeMethods.R | 157 ------------ rstan3/R/rstan.R | 182 -------------- rstan3/R/setAs.R | 81 ------ rstan3/R/stan_csv.R | 432 -------------------------------- rstan3/R/zzz.R | 22 -- 10 files changed, 1478 deletions(-) delete mode 100644 rstan3/.Rbuildignore delete mode 100644 rstan3/DESCRIPTION delete mode 100644 rstan3/NAMESPACE delete mode 100644 rstan3/R/AllClass.R delete mode 100644 rstan3/R/StanFitMCMCMethods.R delete mode 100644 rstan3/R/StanTypeMethods.R delete mode 100644 rstan3/R/rstan.R delete mode 100644 rstan3/R/setAs.R delete mode 100644 rstan3/R/stan_csv.R delete mode 100644 rstan3/R/zzz.R diff --git a/rstan3/.Rbuildignore b/rstan3/.Rbuildignore deleted file mode 100644 index f4bfba209..000000000 --- a/rstan3/.Rbuildignore +++ /dev/null @@ -1,23 +0,0 @@ -inst/unitTests/.*\.html$ -inst/unitTests/.*\.txt$ -vignettes/makefile$ -vignettes/.*\.tex$ -_vignettes/.*\.pdf$ -vignettes/.*\.aux$ -vignettes/.*\.log$ -vignettes/.*\.out$ -vignettes/.*\.bbl$ -vignettes/.*\.toc$ -tests/.*[.]html$ -tests/.*[.]txt$ -\.gitignore$ -Rbuildignore.*$ -.*\.new$ -.*\.old$ -.*\.tmp$ -.*\.sh$ -.*\.so$ -.*\.o$ -.*\.out$ -.*\.backup$ - diff --git a/rstan3/DESCRIPTION b/rstan3/DESCRIPTION deleted file mode 100644 index 25ec08a14..000000000 --- a/rstan3/DESCRIPTION +++ /dev/null @@ -1,26 +0,0 @@ -Package: rstan3 -Type: Package -Title: A Provisional R Interface to Stan -Version: 2.16.0 -Date: 2017-05-22 -Authors@R: c(person("Jiqiang", "Guo", email = "guojq28@gmail.com", role = "aut"), - person("Daniel", "Lee", role = "ctb"), - person("Ben", "Goodrich", email = "benjamin.goodrich@columbia.edu", role = c("cre", "aut"))) -Description: This is a preview of version 3.x of the rstan package. Do not use it. -License: GPL (>=3) -Imports: - Rcpp (>= 0.12.0), - rstan (>= 2.16.2), - sourcetools, - stats, - stats4, - utils -Depends: - R (>= 3.0.2), - methods -Enhances: - rstan -Suggests: - roxygen2 -URL: http://mc-stan.org -RoxygenNote: 5.0.1 diff --git a/rstan3/NAMESPACE b/rstan3/NAMESPACE deleted file mode 100644 index 7b0d9b684..000000000 --- a/rstan3/NAMESPACE +++ /dev/null @@ -1,12 +0,0 @@ -# Generated by roxygen2: do not edit by hand - -export(stan_build) -importFrom(Rcpp,".__C__C++Object") -importFrom(methods,as) -importFrom(methods,callGeneric) -importFrom(methods,new) -importFrom(methods,setRefClass) -importFrom(stats,fivenum) -importFrom(utils,head) -importFrom(utils,str) -importFrom(utils,tail) diff --git a/rstan3/R/AllClass.R b/rstan3/R/AllClass.R deleted file mode 100644 index 85529d899..000000000 --- a/rstan3/R/AllClass.R +++ /dev/null @@ -1,450 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2015, 2016, 2017 Trustees of Columbia University -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -help_from_instance <- function(topic) { - "Returns documentation string for the $method given by 'topic'" - CLASS <- .self$getRefClass() - if (missing(topic)) print(CLASS$help()) - else { - if (is.name(substitute(topic))) topic <- substitute(topic) - print(CLASS$help(parse(text = topic))) - } - return(invisible(NULL)) -} - -#' A ReferenceClass representing a Stan module -#' @importFrom methods as callGeneric new setRefClass -#' @importFrom Rcpp .__C__C++Object -#' @importFrom stats fivenum -#' @importFrom utils head str tail -StanModule <- setRefClass("StanModule", fields = list(module = "C++Object"), - methods = list( -help = help_from_instance, -hmc_nuts_diag_e_adapt = function(chains = 4, iter = 2000, warmup = 1000, - thin = 1, init = "random", - sample_file = NULL, diagnostic_file = NULL, - adapt_gamma = 0.05, adapt_delta = 0.8, - adapt_kappa = 0.75, adapt_t0 = 10, - adapt_init_buffer = 75, - adapt_term_buffer = 50, - adapt_window = 25, stepsize = 1, - stepsize_jitter = 0, max_treedepth = 10, - ...) { - "No U-Turn Sampling (NUTS) with adaptive integration time and a Euclidean metric - \n'chains' = 4: number of Markov chains to execute - \n'iter' = 2000: total number of iterations to execute per chain, including 'warmup' - \n'warmup' = 1000: number of warmup iterations - \n'thin' = 1: thinning interval - \n'sample_file' = NULL: optional path specifying where to store draws on disk - \n'diagnostic_file' = NULL: optional path specifying where to store diagnostics on disk - \n'adapt_gamma' = 0.05: gamma tuning parameter - \n'adapt_delta' = 0.8: delta tuning parameter - \n'adapt_kappa' = 0.75: kappa tuning parameter - \n'adapt_t0' = 10: t0 tuning parameter - \n'adapt_init_buffer' = 75: pre-adaptation buffer - \n'adapt_term_buffer' = 50: final adaptation buffer - \n'adapt_window' = 25: interior phase buffer - \n'stepsize' = 1: initial stepsize, final stepsize is adapted - \n'stepsize_jitter' = 0: bad idea, do not change this - \n'max_treedepth' = 10: maximum number of leapfrog steps per iteration - \n'...' = various other arguments such as: - \nchain_id (integer >= 0) - \ninit_r (numeric > 0) - \nrefresh (integer >= 0) - \nsave_warmup (logical)" - - control <- list(adapt_engaged = TRUE, adapt_gamma = adapt_gamma, - adapt_kappa = adapt_kappa, adapt_t0 = adapt_t0, - adapt_init_buffer = as.integer(adapt_init_buffer), - adapt_term_buffer = as.integer(adapt_term_buffer), - adapt_window = as.integer(adapt_window), - stepsize = stepsize, stepsize_jitter = stepsize_jitter, - max_treedepth = as.integer(max_treedepth)) - args_list <- rstan:::config_argss(chains = as.integer(chains), - iter = as.integer(iter), - warmup = as.integer(warmup), - thin = as.integer(thin), init, - seed = 12345L, # FIXME - sample_file, diagnostic_file, - algorithm = "NUTS", - control = control, ...) - samples <- lapply(args_list, FUN = function(x) .self$module$call_sampler(x)) - return(StanFitMCMC(sample_draws = samples)) -}, -sample = function() { - "calls the 'hmc_nuts_diag_e_adapt' method with its default arguments" - return(.self$hmc_nuts_diag_e_adapt()) -} -) -) - -# A ReferenceClass representing the fixed aspects of a Stan program -setRefClass("StanBuild", - contains = "VIRTUAL", - fields = list(program.name = "character", - shared.object = "character", - seed = "integer"), - # in child classes there is one field for the C++ class and - # fields for each element of the data block in a Stan program - methods = list( -finalize = function() { - try(dyn.unload(.self$shared.object), silent = TRUE) - return(invisible(NULL)) -}, -# FIXME: Should this be an initialize method? -data = function(data, search.env = parent.env(.self), inherits.env = FALSE) { - fields <- tail(names(.self$getRefClass()$fields()), -4L) - if (length(data) == 0) { - data <- mget(fields, envir = search.env, inherits = inherits.env, - ifnotfound = list(function(...) NULL)) - } - nms <- names(data) - for (i in seq_along(data)) if (!is.null(data[[i]])) { - if (nms[i] %in% fields) - .self$field(nms[i], as(data[[i]], class(.self$field(nms[i])))) - else - warning(nms[i], " is not declared in the data block of the Stan program") - } -}, -# Question: Should build be a method or a standalone function() -build = function() { - nms <- names(.self$getRefClass()$fields())[-c(1:3)] - dummy <- function(...) stop("This function should not get called") - datlist <- sapply(nms[-1], simplify = FALSE, - FUN = function(nm) .self$field(nm)@values) - mod <- new(.self$field(nms[1]), datlist, .self$seed, dummy) - # FIXME: why the hell does this require rstan3::: ? - return(rstan3:::StanModule(module = mod)) -}, -show = function() { - nms <- tail(names(.self$getRefClass()$fields()), -3L) - cat("Stan program", sQuote(.self$program.name), "with data", "\n", sep = " ") - for (nm in nms) { - cat(nm, ": ", sep = "") - print(.self$field(nm)) - cat("\n") - } -}, -help = help_from_instance -)) - -#' Constructor -#' -#' @param ... Objects passed to the data block of the Stan program -#' @param file A character vector of length one specifying a path to a file -#' written in the Stan language, which is ignored if \code{model_code} is -#' specified -#' @param model_name A character vector of length one giving a name to the -#' Stan program, which is autogenerated by default -#' @param model_code An optional character vector written in the Stan language -#' @param includes An optional character vector specifying includes and stuff -#' @param seed An integer for Stan's PRNG seed -#' @export -stan_build <- function(..., - file = file.choose(), - model_name = NA_character_, - model_code, - includes = character(), - seed = sample.int(.Machine$integer.max, 1L)) { - allow_undefined <- length(includes) > 0 - if (missing(model_code)) { - md5 <- tools::md5sum(file) - tf <- file.path(tempdir(), paste0(md5, ".stan")) - file.copy(from = file, to = tf) - if (is.na(model_name)) - model_name <- sub("\\..*$", "", basename(file)) - } - else { - tf <- tempfile(fileext = ".stan") - writeLines(model_code, tf) - md5 <- tools::md5sum(tf) - new_name <- file.path(tempdir(), paste0(md5, ".stan")) - file.rename(from = tf, to = new_name) - tf <- new_name - if (is.na(model_name)) - model_name <- "anon_model" - } - on.exit(unlink(tf)) - ret <- rstan::stanc_builder(tf, verbose = FALSE, - allow_undefined = allow_undefined) - class_name <- paste("model", ret$model_name, sep = "_") - inc <- paste("#define STAN__SERVICES__COMMAND_HPP\n", - "// [[Rcpp::depends(rstan)]]\n", - if(length(includes)) ret$cppcode else - sub("(class[[:space:]]+[A-Za-z_][A-Za-z0-9_]*[[:space:]]*: public prob_grad \\{)", - paste(includes, "\\1"), ret$cppcode), - "#include \n", - rstan:::get_Rcpp_module_def_code(class_name), - sep = '') - stuff <- Rcpp::sourceCpp(code = inc, env = environment()) - - cppcode <- scan(what = character(), sep = "\n", quiet = TRUE, - text = ret$cppcode) - private <- grep("^private:$", cppcode) + 1L - public <- grep("^public:$", cppcode) - 1L - objects <- gsub("^.* ([0-9A-Za-z_]+).*;.*$", "\\1", cppcode[private:public]) - in_data <- grep("context__.vals_", cppcode, fixed = TRUE, - value = TRUE) - in_data <- gsub("^.*\\(\"(.*)\"\\).*;$", "\\1", in_data) - objects <- intersect(objects, in_data) - - fields <- make_fields(ret$model_code, objects) - fields <- c("C++Class", fields) - obj_name <- paste0("stan_fit4", class_name) - names(fields)[1] <- obj_name - SPWD <- setRefClass(md5, contains = "StanBuild", fields = fields, - where = .GlobalEnv) - SPWD$lock(c("shared.object", obj_name, "seed")) - out <- SPWD(program.name = model_name, - shared.object = dir(stuff$buildDirectory, pattern = "\\.so$", - full.names = TRUE), - seed = as.integer(seed)) - out$field(obj_name, get(obj_name, inherits = FALSE)) - m <- match.call() - dotlist <- list(...) - no_names <- is.null(names(dotlist)) - has_name <- if (no_names) FALSE else nzchar(names(dotlist)) - m$file <- m$model_name <- m$model_code <- m$includes <- m$seed <- NULL - nms <- as.character(m)[-1L] - if (no_names) names(dotlist) <- nms - else names(dotlist)[!has_name] <- nms[!has_name] - out$data(dotlist) - return(out$build()) -} - -# Dynamically creates field names corresponding to objects in the data block -make_fields <- function(stan_code, objects) { - if (length(objects) == 0) return(list()) - tokens <- sourcetools::tokenize_string(stan_code) - matches <- sapply(objects, FUN = function(x) which(tokens$value == x)[1]) - rows <- tokens$row[matches] - tokens <- tokens[tokens$row %in% rows,] - semis <- which(tokens$type == "semi") - arrays <- sapply(semis, FUN = function(s) { - while(TRUE) { - s <- s - 1L - if (tokens[s, "value"] == "]") return(TRUE) - if (tokens[s, "type"] == "symbol") return(FALSE) - } - }) - if (length(arrays) == 0) arrays <- integer() - StanTypes <- names(getClass("StanType")@subclasses) - types <- sapply(tokens$value, FUN = match, table = StanTypes) - types <- types[!is.na(types)] - nms <- sapply(1:length(semis), FUN = function(i) { - s <- semis[i] - if (!is.na(arrays[i]) && arrays[i]) while(TRUE) { - s <- s - 1L - if (tokens[s, "value"] == "[") return(tokens[s - 1L, "value"]) - } - else while(TRUE) { - s <- s - 1L - if (tokens[s, "type"] == "symbol") return(tokens[s, "value"]) - } - }) - names(types) <- nms - return(sapply(types, simplify = FALSE, FUN = function(j) StanTypes[j])) -} - -setClassUnion("scalarORarray", c("numeric", "array")) - -#' A S4 class tree mirroring the types in the Stan language -#' -#' These classes are not intended to be generated by users, but there are -#' associated S4 methods that can be called by users. -#' -#' The \code{StanType-class} is virtual and the inheritance logic is as -#' follows: -#' \itemize{ -#' \item real: inherits from \code{StanType-class} and represents -#' real scalars -#' \itemize{ -#' \item int: inherits from \code{real-class} and represents -#' integer scalars with optional slots for 'levels' and 'ordered' -#' that can be used when coercing a \R \code{\link{factor}} -#' } -#' \item mat: inherits from \code{StanType-class} and represents -#' real matrices -#' \itemize{ -#' \item cov_matrix: inherits from \code{mat-class} and -#' represents covariance matrices -#' \itemize{ -#' \item corr_matrix: inherits from \code{mat-class} and -#' represents correlation matrices -#' } -#' \item cholesky_factor_cov: inherits from \code{mat-class} and -#' represents a Cholesky factor of a covariance matrix -#' \itemize{ -#' \item cholesky_factor_corr: inherits from -#' \code{cholesky_factor_cov-class} and represents a Cholesky -#' factor of a correlation matrix -#' } -#' \item vec: inherits from \code{mat-class} and represents -#' a vector, even though in the Stan language a vector and a -#' matrix are both primitive types -#' \itemize{ -#' \item row_vector: inherits from \code{vec-class} and -#' represents a row vector, even though in the Stan language a -#' vector and a row_vector are both primitive types -#' \item simplex: inherits from \code{vec-class} and -#' represents a simplex -#' \item unit_vector: inherits from \code{vec-class} and represents -#' a unit vector on a hypersphere -#' \item ordered: inherits from \code{vec-class} and represents -#' an ordered vector -#' \itemize{ -#' \item positive_ordered: inherits from \code{ordered-class} and -#' represents an ordered vector whose elements are positive -#' } -#' } -#' } -#' } -#' -#' @slot values A numeric (possibly integer) scalar or \code{link[base]{array}} -#' @slot lower A numeric scalar lower bound -#' @slot upper A numeric scalar upper bound -#' @slot type A length-one character vector among -#' \itemize{ -#' \item data -#' \item parameter -#' \item transformed parameter -#' \item generated quantity -#' \item diagnostic -#' } -#' @slot constrained A logical scalar indicating whether an unknown is in -#' its constrained state -#' -setClass("StanType", contains = "VIRTUAL", - slots = list(values = "scalarORarray", - lower = "numeric", upper = "numeric", - type = "character", constrained = "logical"), - prototype = list(values = array(NA_integer_, dim = 0L), - lower = -Inf, upper = Inf, - type = "data", - constrained = TRUE), - validity = function(object) { - if (!is.numeric(object@values)) - return("'values' must be numeric") - if (length(object@type) != 1) - return("'type' must be of length one") - if (length(object@lower) != 1) - return("'lower' must be of length one") - if (length(object@upper) != 1) - return("'upper' must be of length one") - if (object@lower > object@upper) - return("'lower' must be less than or equal to upper") - if (any(object@values < object@lower)) - return("'lower' bound violated") - if (any(object@values > object@upper)) - return("'upper' bound violated") - types <- c("data", "parameter", "transformed parameter", - "generated quantity", "diagnostic") - if (!(object@type %in% types)) - return(paste("'type' must be among:", - paste(types, collapse = "\n"), sep = "\n")) - return(TRUE) - } -) - -#' @rdname StanType-class -setClass("real", contains = "StanType") - -#' @rdname StanType-class -setClass("int", contains = "real", - slots = list(labels = "character", ordered = "logical"), - prototype = list(lower = -.Machine$integer.max, - upper = .Machine$integer.max), - validity = function(object) { - if (!is.integer(object@values)) return("'values' must contain integers") - if (identical(object@type, "paramter") || - identical(object@type, "transformed parameters")) - return("'int' cannot be used in parameters or transformed parameters") - return(TRUE) - } -) - -#' @rdname StanType-class -setClass("mat", contains = "StanType") - -#' @rdname StanType-class -setClass("cov_matrix", contains = "mat") -#' @rdname StanType-class -setClass("corr_matrix", contains = "cov_matrix") -#' @rdname StanType-class -setClass("cholesky_factor_cov", contains = "matrix") -#' @rdname StanType-class -setClass("cholesky_factor_corr", contains = "cholesky_factor_cov") -#' @rdname StanType-class -setClass("vec", contains = "mat") - -#' @rdname StanType-class -setClass("row_vector", contains = "vec") -#' @rdname StanType-class -setClass("simplex", contains = "vec", validity = function(object) { - if (any(object@values < 0)) - return("'values' must be non-negative") - if (!(sum(object@values) %in% 0:1)) - return("'values' must sum to 1") - return(TRUE) - } -) -#' @rdname StanType-class -setClass("unit_vector", contains = "vec", validity = function(object) { - if (any(object@values < -1)) - return("'values' must >= -1") - if (any(object@values > 1)) - return("'values' must <= 1") - if (!(sum(object@values ^ 2) %in% 0:1)) - return("sum of squared 'values' must be 1") - return(TRUE) - } -) -#' @rdname StanType-class -setClass("ordered", contains = "vec", validity = function(object) { - if (any(diff(object@values) < 0)) - return("'values' must be in increasing order") - return(TRUE) - } -) -#' @rdname StanType-class -setClass("positive_ordered", contains = "ordered", validity = function(object) { - if (any(object@values < 0)) return("'values' must be non-negative") - return(TRUE) - } -) - -not_implemented <- function() stop("this is not implemented yet") - -#' A ReferenceClass representing MCMC output -#' -#' @field warmup_draws A list of \code{\link{StanType}}s during warmup -#' @field sample_draws A list of \code{\link{StanType}}s after warmup -#' -#' @seealso \code{\link{StanFitOptimization-methods}} -StanFitMCMC <- - setRefClass("StanFitMCMC", fields = list(warmup_draws = "list", - sample_draws = "list"), - # FIXME: add other extraction methods but without the get_ prefixes - methods = list(show = not_implemented, - summary = not_implemented, - as.mcmc.list = not_implemented, - extract = not_implemented, - add_params = not_implemented, - help = help_from_instance) -) -StanFitMCMC$lock(setdiff(names(StanFitMCMC$fields()), "added_draws")) - diff --git a/rstan3/R/StanFitMCMCMethods.R b/rstan3/R/StanFitMCMCMethods.R deleted file mode 100644 index dc9ca33bd..000000000 --- a/rstan3/R/StanFitMCMCMethods.R +++ /dev/null @@ -1,93 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2015 Jiqiang Guo and Benjamin Goodrich -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -StanFitMCMC$methods(show = function() { - "Show a brief version of the posterior" - print("brief information about the fit") - # dims <- dim(.self) - # cat(dims[1], "unknowns,", dims[2], "chains, and", - # dims[3], "retained samples") - # return(invisible(NULL)) -}) - -StanFitMCMC$methods(summary = function() { - "Compute the summary of the posterior" - print("detailed information about the fit") - # out <- sapply(.self$sample_draws, FUN = stats4::summary, - # simplify = FALSE) - # return(t(do.call(cbind, args = out))) -}) - -StanFitMCMC$methods(as.mcmc.list = function() { - "Convert to an mcmc.list object for use with the coda package" - arr <- as.array(.self) - return(coda::as.mcmc.list(lapply(1:ncol(arr), FUN = function(j) { - mcmc(t(arr[,j,])) - }))) -}) - -# S3 methods for StanFitMCMC -"[.StanFitMCMC" <- function(x, i, ...) { - if (length(i) == 1) { - out <- x$sample_draws[[i]] - } - else { - out <- x$sample_draws[i] - } - return(out) -} - -"[[.StanFitMCMC" <- function(x, i, ...) { - if (is.character(i)) { - return(sapply(i, simplify = length(i) == 1, FUN = function(y) { - param <- grep(paste0("^", y, "["), names(x@sample_params)) - return(x@sample_params[[param]][i,]) - })) - } - else if (is.numeric(i)) { - breaks <- cumsum(sapply(x@sample_params, nrow)) - return(sapply(i, simplify = length(i) == 1, FUN = function(y) { - param <- which(i <= breaks)[1] - return(x@sample_params[[param]][i - param + 1,]) - })) - } -} - -dim.StanFitMCMC <- function(x) { - dims <- sapply(x$sample_draws, FUN = dim, simplify = FALSE) - params <- sum(unlist(sapply(dims, FUN = head, n = -2L))) - chains <- tail(dims[[1]], 2) - iterations <- chains[2] - chains <- chains[1] - return(c(params, chains, iterations)) -} - -as.array.StanFitMCMC <- function(x) { - out <- do.call(rbind, args = - sapply(x$sample_draws, FUN = as.matrix, simplify = FALSE)) - dim(out) <- c(nrow(out), dim(x)[-1]) - return(out) -} - -as.matrix.StanFitMCMC <- function(x) { - out <- as.array(x) - n <- nrow(out) - dim(out) <- c(n, length(out) %/% n) - return(out) -} - -names.StanFitMCMC <- function(x) names(x$sample_draws) diff --git a/rstan3/R/StanTypeMethods.R b/rstan3/R/StanTypeMethods.R deleted file mode 100644 index 03caaa8ff..000000000 --- a/rstan3/R/StanTypeMethods.R +++ /dev/null @@ -1,157 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2015 Jiqiang Guo and Benjamin Goodrich -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -setMethod("show", signature = "StanType", definition = function(object) { - if (length(object@values)) cat(str(object@values)) - else cat(class(object)) -}) - -setMethod("summary", signature = "StanType", definition = function(object) { - leading <- head(dim(object), n = -2L) - out <- apply(object@values, MARGIN = leading, FUN = fivenum) - dim(out) <- c(5L, length(out) %/% 5L) - return(out) -}) - -# setMethod("show", signature = "StanFactor", definition = function(object) { -# if( length(dim(object@values)) == 1) { -# cat("StanType is not filled in yet") -# return(invisible(NULL)) -# } -# cat(object@name, "has", tail(dim(object@values), 1), "total draws") -# return(invisible(NULL)) -# }) -# -# setMethod("summary", signature = "StanFactor", definition = function(object) { -# if( length(dim(object@values)) == 1) { -# cat("StanType is not filled in yet") -# return(invisible(NULL)) -# } -# values <- factor(c(object@values), levels = 1:length(object@labels), -# levels = levels, ordered = object@ordered) -# tab <- table(values) -# return(tab) -# }) - -setMethod("summary", signature = "cov_matrix", definition = function(object) { - # do not summarize the upper triangle - dims <- dim(object@values) - len <- length(dims) - rows <- dims[len - 2L] - cols <- dims[len - 1L] - iterations <- dims[len] - chains <- object@chains - iterations_per_chain <- iterations %/% chains - monitor_covariance <- function(r, c = r, x) { - x <- aperm(x, perm = c(iterations, dims[-len])) - dim(x) <- c(iterations_per_chain, chains, dim(x)[-1L]) - stop("FIXME: Implement monitor") - monitor(x[,,r,c]) - } - if (len == 3) { - out <- rbind( sapply(1:rows, FUN = monitor_covariance, x = object@values), - if(rows > 1) sapply(2:rows, FUN = function(r) { - sapply(1:(r-1L), monitor_covariance, r = r, x = object@values) - }) ) - } - else { - out <- apply(object@values, MARGIN = 1:(len-3L), FUN = function(y) { - sapply(1:rows, FUN = monitor_covariance, x = y) - }) - if(rows > 1) out <- rbind(out, apply(object@values, MARGIN = 1:(len-3L), - FUN = function(y) { - sapply(2:rows, FUN = function(r) { - sapply(1:(r - 1L), monitor_covariance(r = r, - x = y)) - }) - })) - } - return(t(out)) -}) - -# (Re)Defining these S4groupGeneric functions facilitates functional programming -setMethod("Ops", signature = signature("ANY", "StanType"), - definition = function(e1, e2) { - return(new(class(e2), values = callGeneric(e1, e2@values))) -}) -setMethod("Ops", signature = signature("StanType", "ANY"), - definition = function(e1, e2) { - return(new(class(e1), values = callGeneric(e1@values, e2))) -}) -setMethod("Ops", signature = signature("StanType", "StanType"), - definition = function(e1, e2) { - return(new(class(e1), values = callGeneric(e1@values, e2@values))) -}) -setMethod("Math", signature = "StanType", - definition = function(x) { - return(new(class(x), values = callGeneric(x@values))) -}) -setMethod("log", signature = "StanType", - definition = function(x, ...) { - return(new(class(x), values = log(x@values, ...))) -}) -setMethod("trunc", signature = "StanType", - definition = function(x, ...) { - return(new(class(x), values = trunc(x@values, ...))) -}) -setMethod("Math2", signature = "StanType", - definition = function(x, digits) { - return(new(class(x), arr = callGeneric(x@values, digits))) -}) -setMethod("Summary", signature = "StanType", - definition = function(x, ..., na.rm = FALSE) { - return(new(class(x), values = callGeneric(x@values, ..., na.rm))) -}) -setMethod("Complex", signature = "StanType", - definition = function(z) { - return(new(class(z), values = callGeneric(z@values))) -}) - -# S3 Methods for a StanType -"[.StanType" <- function(x, i, j, ...) { - if (missing(j)) return(new(class(x), name = rownames(x@values)[i], - values = x@values[i,], type = x@type)) - else return(values = x@values[i,j,...]) -} - -"[.real" <- function(x, i = 1L, j, ...) { - if (is.numeric(i) && i != 1) stop("'i' can only be 1") - if (is.character(i) && i != x@name) - stop(paste("'i' can only be", x@name)) - - if (missing(j)) return(x) - else return(values = x@values[i,j,...]) -} - -dim.StanType <- function(x) dim(x@values) - -as.array.StanType <- function(x, ...) x@values -as.matrix.StanType <- function(x, ...) { - dims <- dim(x) - leading <- head(dims, n = -2L) - rows <- prod(leading) - cols <- prod(tail(dims, n = 2L)) - mat <- x@values - dim(mat) <- c(rows, cols) - combos <- do.call(expand.grid, args = - lapply(leading, function(z) 1:z)) - rownames(mat) <- paste0(x@name, "[", - apply(combos, 1, paste, collapse = ","), "]") - return(mat) -} -as.matrix.real <- function(x, ...) x@values - diff --git a/rstan3/R/rstan.R b/rstan3/R/rstan.R deleted file mode 100644 index 64a1b1dc7..000000000 --- a/rstan3/R/rstan.R +++ /dev/null @@ -1,182 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2015, 2016, 2017 Trustees of Columbia University -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -#' Estimate and diagnose models using Stan -#' -#' This R package provides an interface to the Stan C++ library or more -#' specifically links to the \pkg{StanHeaders} package, which provides the -#' header files of the Stan C++ libraries. -#' -#' Starting from \pkg{rstan} 3.0, the syntax of \code{\link[methods]{ReferenceClasses}} -#' is heavily utilized by the \pkg{rstan} package. This syntax may be unfamiliar -#' to many \R users but is similar to that in Python. In particular, various -#' \code{\link{ReferenceClasses}} are defined with fields and methods. These -#' methods are \code{\link{function}}s that may utilize the fields in addition -#' to any arguments of the methods. Perhaps the biggest syntactic difference is -#' that these methods are called using the \code{$} operator, such as -#' \code{object$plot()} rather than \code{plot(object)}. To list the methods -#' that are defined for \code{\link[methods]{ReferenceClasses}}, postfix \code{$} -#' to the object name and press the Tab key. To obtain help for any defined -#' method, use the \code{$help(topic)} method, where \code{topic} is the name -#' of the method you need help with. -#' -#' The basic steps for defining, estimating, and diagnosing a Stan model are as -#' follows. -#' -#' @section Step 0 --- Write a Stan program: -#' Preferably save a file with a .stan extension to your hard disk that -#' constitutes a valid Stan program or you can write a valid -#' Stan program as a \R \code{\link{character}} vector. A simple illustration -#' of a Stan program is given in the examples section below. More details about -#' the Stan language are at \url{http://mc-stan.org/documentation}, but in brief -#' a Stan program has one or more of the following blocks contained within curly -#' braces: -#' \itemize{ -#' \item functions: an optional section for users to define functions that -#' can subsequently be used in that Stan program only or can be exported -#' to the \R by calling \code{link[rstan]{expose_stan_functions}} -#' \item data: an all-but-required section for users to declare known -#' quantities, including the outcome to be modeled, its predictors, -#' hyperparameters, dimensions, etc. These are passed from R. -#' \item transformed data: an optional section for users to define other -#' known quantities that are related to the objects declared in the \code{data} -#' block and hence is only executed once. It is also possible to do random -#' number generation in this block. -#' \item parameters: an all-but-required section to declare the unknowns -#' that are being estimated and hence are saved in the output -#' \item transformed parameters: an optional section for users to define -#' unknown quantities that are conditionally deterministic given the -#' objects in \code{parameters}, \code{transformed data}, and \code{data}. -#' The \code{transformed parameters} block is executed every time the -#' objective function is evaluated, the constraints on the objects therein -#' are checked, and the objects are saved in the output -#' \item model: an all-but-required section for users to build up an -#' expression for the objective function, which may consist of both priors -#' and a likelihood to form a posterior kernel. Local objects can be declared -#' at the top of the \code{model} block but their are not saved in the output. -#' \item generated quantities: an optional section where users can define -#' additional objects that depend on anything in the previous sections -#' \emph{except} the \code{model} section and may involve random-number -#' generation. The \code{generated quantities} section is executed once -#' per iteration and is useful for producing posterior predictive -#' distributions. -#' } -#' The Stan language supports the several types of objects. These same types are now -#' available in \R via S4 classes, as documented by the \code{\link{StanType-class}}. -#' -#' @section Step 1 --- Build your Stan program: -#' Either pass a file path to the \code{file} argument of the -#' \code{\link{stan_config}} constructor or pass specify the \code{mode_code} -#' argument. For example, \code{post <- stan_build(file = "MyStanFile.stan")} -#' would locate a file called MyStanFile.stan in the working directory on the -#' hard disk. If the \code{file} argument is unspecified, a dialogue box will -#' open (in \code{\link[base]{interactive}} mode) that will prompt you to a file -#' on your hard disk. An illustration of the \code{model_code} argument is given -#' in the examples section below. Either approach will cause the Stan syntax to -#' be parsed to C++ and (if successful) compiled, which may take a considerable -#' amount of time. This will create an object of \code{\link{StanModule-class}}. -#' -#' By default \code{stan_build} will look in the calling environment for R -#' objects named in the data block of the Stan program and will use them to -#' instantiate the C++ object that ultimately estimates the unknown parameters. -#' However, you can (and should) specify these arguments explicitly when calling -#' \code{stan_build} as illustrated in the example below. If the R objects have -#' the same names as the Stan objects, then you can simply pass them. Otherwise, -#' use the \code{tag = object} R syntax as for any other \code{\link{list}}, -#' where \code{tag} is the name of the object in the data block of the Stan -#' program. -#' -#' @section Step 2 --- Estimate the unknowns your model: -#' There are a variety of methods for an object of \code{\link{StanModule-class}} -#' that can be invoked to estimate the unknowns of your model in various ways. -#' The preferred way is to draw from the posterior distribution using Markov -#' Chain Monte Carlo (MCMC) sampling, but it is also possible to use Automatic -#' Differentiation Variational Inference (ADVI) with either a meanfield or a -#' fullrank assumption on the covariance matrix of the multivariate normal -#' distribution that is closest to the posterior distribution in the -#' unconstrained space. Finally, it is possible to find the mode of a -#' log-likelihood (if not priors are specified) or a posterior mode (if priors -#' are specified) using various deterministic optimization algorithms that also -#' make use of automatic differentiation to evaluate the gradient with respect -#' to the unknowns in the unconstrained parameter space. These algorithms can -#' be called with their default arguments by invoking one of the following -#' methods: -#' \itemize{ -#' \item \code{sample}: MCMC, specifically No U-Turn Sampling (NUTS) with -#' adaptive integration time, a tuned diagonal mass matrix, and a -#' Euclidean metric -#' \item \code{advi}: ADVI with a diagonal covariance matrix -#' \item \code{maximize}: The LBFGS mode-finding algorithm -#' } -#' Other algorithms can be invoked (possibly with non-default arguments) via -#' the corresponding methods, which are documented with the -#' \code{\link{StanModule-class}}. -#' -#' @section Step 3 --- Diagnose any problems: -#' The best way to diagnose problems with MCMC is to call the -#' \code{\link[shinyStan]{launch_shinystan}} function in the -#' \pkg{shinyStan} package, but this requires user interactivity. Many -#' static plots can be produced interactively or by a script (via the -#' \pkg{bayesplot} package) by invoking the methods for an object of -#' \code{\link{StanMCMCFit-class}}. However, the most basic methods are -#' \code{show} and, for more detail, \code{summary}. -#' -#' @examples -#' # create data in R -#' J <- 8 # number of schools -#' y <- c(28, 8, -3, 7, -1, 1, 18, 12) # estimated treatment effects -#' sigma <- c(15, 10, 16, 11, 9, 11, 10, 18) # estimated standard errors -#' -#' # Step 0 --- Write a Stan program -#' mc <- ' -#' data { /* objects conditioned on in Bayes Law */ -#' int J; // syntax for an non-negative integer scalar -#' real y[J]; // syntax for an array of real numbers -#' vector[J] sigma; // syntax for a non-negative vector of reals -#' } -#' parameters { /* variables whose posterior distribution is sought */ -#' real mu; -#' real tau; -#' vector[J] eta; -#' } -#' transformed parameters { /* conditionally deterministic variables */ -#' vector[J] theta = mu + tau * eta; -#' } -#' model { /* build up expression for the posterior distribution */ -#' // implicit: improper priors on mu and tau -#' target += normal_lpdf(eta | 0, 1); -#' // normal prior on eta implies: theta ~ normal(mu, tau) -#' target += normal_lpdf(y | theta, sigma); // normal likelihood -#' } -#' ' -#' \donttest{ -#' # Step 1 --- Build your Stan program -#' # a path to file is preferable but a character vector works -#' post <- stan_build(model_code = mc, model_name = "8schools", seed = 12345, -#' data = list(J = J, y = y, sigma = sigma)) -#' -#' # Step 2 --- Estimate the unknowns ofyour model -#' post$help("hmc_nuts_diag_e_adapt") # this is what is called by $sample() -#' draws <- post$sample() # draws from posterior distribution -#' -#' # Step 3 --- Diagnose any problems -#' draws # shows basic characteristics of the posterior distribution -#' draws$summary() # shows more details about the posterior distribution -#' } -#' @docType package -#' @name rstan3 -NULL diff --git a/rstan3/R/setAs.R b/rstan3/R/setAs.R deleted file mode 100644 index b011dd74d..000000000 --- a/rstan3/R/setAs.R +++ /dev/null @@ -1,81 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2017 Trustees of Columbia University -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -setAs("numeric", "real", function(from) { - new("real", values = from) -}) -setAs("list", "real", function(from) { - from <- simplify2array(from) - dims <- dim(from) - len <- length(dims) - if (len > 1L) - from <- aperm(from, perm = c(dims[len], dims[-len])) - as(from, "real") -}) -setAs("matrix", "real", function(from) { - new("real", values = as.array(from)) -}) - -setAs("factor", "int", function(from) { - new("int", values = as.integer(from), - labels = levels(from), ordered = is.ordered(from)) -}) -setAs("numeric", "int", function(from) { - new("int", values = as.integer(from)) -}) -setAs("integer", "int", function(from) { - new("int", values = from) -}) -setAs("array", "int", function(from) { - new("int", values = from) -}) - -setAs("list", "int", function(from) { - from <- simplify2array(from) - dims <- dim(from) - len <- length(dims) - if (len > 1L) - from <- aperm(from, perm = c(dims[len], dims[-len])) - as(from, "int") -}) - -setAs("matrix", "mat", function(from) { - new("mat", values = from) -}) -setAs("list", "mat", function(from) { - from <- simplify2array(from) - dims <- dim(from) - len <- length(dims) - # FIXME matrices with 1 row or column - from <- aperm(from, perm = c(dims[len], dims[-len])) - as(from, "mat") -}) - -setAs("numeric", "vec", function(from) { - dims <- dim(from) - if (is.null(dims)) dims <- length(from) - from <- array(from, dim = dims) - new("vec", values = from) -}) -setAs("list", "vec", function(from) { - from <- simplify2array(from) - dims <- dim(from) - len <- length(dims) - if (len > 1L) - from <- aperm(from, perm = c(dims[len], dims[-len])) - as(from, "vec") -}) diff --git a/rstan3/R/stan_csv.R b/rstan3/R/stan_csv.R deleted file mode 100644 index bca920872..000000000 --- a/rstan3/R/stan_csv.R +++ /dev/null @@ -1,432 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017 Trustees of Columbia University -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -#' Read CSV files created by CmdStan or RStan -#' -#' @param csvfiles A character vector of paths to CSV files of MCMC output -#' @return An object of \code{\link{StanFitMCMC-class}} -#' @seealso \code{\link{StanFitMCMC-class}} -read_stan_csv <- function(csvfiles) { - - col_major <- TRUE - paridx_fun <- function(names) { - # Args: - # names: names (character vector) such as lp__, treedepth__, stepsize__, - # alpha, beta.1, - # Returns: - # The indexes in the names that are parameters other than lp__, - # treedepth__, or stepsize__. The vector has attribute meta - # with the indexes of 'treedepth__', 'lp__', and 'stepsize__' - # if available. - - sampler_param_names <- c('lp__', 'accept_stat__', 'treedepth__', 'stepsize__', 'divergent__', 'n_leapfrog__') - metaidx <- match(sampler_param_names, names) - names(metaidx) <- sampler_param_names - paridx <- setdiff(seq_along(names), metaidx) - attr(paridx, "meta") <- metaidx[!sapply(metaidx, is.na)] - paridx - } - - read_comments <- function(f, n = -1) { - # Read comments beginning with `#` - # Args: - # f: the filename - # n: max number of line; -1 means all - # Returns: - # a vector of strings - con <- file(f, 'r') - comments <- list() - iter <- 0 - while (length(input <- readLines(con, n = 1)) > 0) { - if (n > 0 && n <= iter) break; - if (grepl("#", input)) { - comments <- c(comments, gsub("^.*#", "#", input)) - iter <- iter + 1 - } - } - close(con) - do.call(c, comments) - } - - read_csv_header <- function(f, comment.char = '#') { - # Read the header of a csv file (the first line not beginning with - # comment.char). And the line number is return as attribute of name 'lineno'. - con <- file(f, 'r') - niter <- 0 - while (length(input <- readLines(con, n = 1)) > 0) { - niter <- niter + 1 - if (!grepl(comment.char, input)) break; - } - header <- input - attr(header, "lineno") <- niter - close(con) - header - } - - parse_stancsv_comments <- function(comments) { - # Parse the comments in Stan CSV files to get information such as - # iter, thin, seed, etc. This is specific to the CSV files - # generated from Stan - - adapt_term_lineno <- which(grepl("Adaptation terminated", comments)) - time_lineno <- which(grepl("Elapsed Time", comments)) - has_time <- length(time_lineno) > 0 - len <- length(comments) - if (length(adapt_term_lineno) < 1) - adapt_term_lineno <- len - if (length(time_lineno) < 1) - warning("line with \"Elapsed Time\" not found") - - if (adapt_term_lineno == len) - adaptation_info <- '' - else { - if (has_time) - adaptation_info <- paste(comments[(adapt_term_lineno+1):(time_lineno-1)], collapse = '\n') - else - adaptation_info <- paste(comments[(adapt_term_lineno+1):len], collapse = '\n') - } - if (has_time) - time_info <- comments[time_lineno:len] - else - time_info <- '' - comments <- comments[1:(adapt_term_lineno - 1)] - - has_eq <- sapply(comments, function(i) grepl('=', i)) - comments <- comments[has_eq] - comments <- gsub('^#+\\s*|\\s*|\\(Default\\)', '', comments) - eq_pos <- regexpr("=", comments, fixed = TRUE) - names0 <- substr(comments, 0, eq_pos - 1) - values <- as.list(substring(comments, eq_pos + 1)) - - id_idx <- which("id" == names0) - if (length(id_idx) > 0) - names0[id_idx] <- "chain_id" - - compute_iter <- FALSE - id_warmup <- which("num_warmup" == names0) - if (length(id_warmup) > 0) { - names0[id_warmup] <- "warmup" - compute_iter <- TRUE - } - - id_numsamples <- which("num_samples" == names0) - if (length(id_numsamples) > 0) { - names0[id_numsamples] <- "iter" - } - names(values) <- names0; - - add_lst <- list(adaptation_info = adaptation_info, - has_time = has_time, - time_info = time_info) - - sampler_t <- NULL - if (!is.null(values$algorithm) && is.null(values$sampler_t)) { - if (values$algorithm == 'rwm' || values$algorithm == 'Metropolis') - sampler_t <- "Metropolis" - else if (values$algorithm == 'hmc') { - if (values$engine == 'static') sampler_t <- "HMC" - else { - if (values$metric == 'unit_e') sampler_t <- "NUTS(unit_e)" - else if (values$metric == 'diag_e') sampler_t <- "NUTS(diag_e)" - else if (values$metric == 'dense_e') sampler_t <- "NUTS(dense_e)" - } - } - add_lst <- c(add_lst, sampler_t = sampler_t) - } - names1 <- intersect(c("thin", "iter", "warmup", "chain_id", "max_depth", - "num_samples", "num_warmup", "id", - "max_treedepth", "save_warmup"), names0) - names2 <- intersect(c("stepsize", "stepsize_jitter", "adapt_gamma", "adapt_kappa", - "adapt_delta", "gamma", "kappa", "delta", "t0", - "adapt_t0"), names0) - for (z in names1) values[[z]] <- as.integer(values[[z]]) - for (z in names2) values[[z]] <- as.numeric(values[[z]]) - if (compute_iter) values[["iter"]] <- values[["iter"]] + values[["warmup"]] - c(values, add_lst) - } - - unique_par <- function(fnames) { - # obtain parameters from flat names in format of say alpha.1, - # alpha.2, beta.1.1, ..., beta.3.4, --- in this case, return - # c('alpha', 'beta') - unique(gsub('\\..*', '', fnames)) - } - - get_dims_from_fnames <- function(fnames, pname) { - # Get the dimension for a parameter from - # the flatnames such as "alpha.1.1", ..., "alpha.3.4", the - # format of names in the CSV files generated by Stan. - # Currently, this function assume fnames are correctly given. - # Args: - # fnames: a character of names for one (vector/array) parameter - # pname: the name for this vector/array parameter such as "alpha" - # for the above example - - if (missing(pname)) pname <- gsub('\\..*', '', fnames[1]) - - if (length(fnames) == 1 && fnames == pname) - return(integer(0)) # a scalar - - idxs <- sub(pname, '', fnames, fixed = TRUE) - lp <- gregexpr('\\d+', idxs) - - tfun <- function(name, start, i) { - last <- attr(start, 'match.length')[i] + start[i] - as.integer(substr(name, start[i], last)) - } - - dim_len <- length(lp[[1]]) - dims <- integer(dim_len) - for (i in 1:dim_len) { - dimi <- mapply(tfun, idxs, lp, MoreArgs = list(i = i), USE.NAMES = FALSE) - dims[i] <- max(dimi) - } - dims - } - - get_time_from_csv <- function(tlines) { - # get the warmup time and sample time from the commented lines - # about time in the CSV files - # Args: - # tlines: character vector of length 3 (or 2 since the last one is not used) - # from the CSV File. For example, it could be - # # Elapsed Time: 0.005308 seconds (Warm-up) - # 0.003964 seconds (Sampling) - # 0.009272 seconds (Total) - t <- rep(NA, 2) - names(t) <- c("warmup", "sample") - if (length(tlines) < 2) return(t) - warmupt <- gsub(".*#\\s*Elapsed.*:\\s*", "", tlines[1]) - warmupt <- gsub("\\s*seconds.*$", "", warmupt) - samplet <- gsub(".*#\\s*", "", tlines[2]) - samplet <- gsub("\\s*seconds.*$", "", samplet) - t[1] <- as.double(warmupt) - t[2] <- as.double(samplet) - t - } - - all_int_eq <- function(is) { - # tell if all integers in 'is' are the same - if (!all(is.integer(is))) - stop("not all are integers") - min(is) == max(is) - } - - if (length(csvfiles) < 1) - stop("csvfiles does not contain any CSV file name") - - g_skip <- 10 - g_max_comm <- -1 # to read all - - cs_lst <- lapply(csvfiles, function(csv) read_comments(csv, n = g_max_comm)) - cs_lst2 <- lapply(cs_lst, parse_stancsv_comments) - - ss_lst <- vector("list", length(csvfiles)) - for (i in seq_along(csvfiles)) { - header <- read_csv_header(csvfiles[i]) - lineno <- attr(header, 'lineno') - vnames <- strsplit(header, ",")[[1]] - m <- matrix(scan(csvfiles[i], skip = lineno, comment.char = '#', sep = ',', quiet = TRUE), - ncol = length(vnames), byrow = TRUE) - ss_lst[[i]] <- as.data.frame(m) - colnames(ss_lst[[i]]) <- vnames - } - - ## read.csv is slow for large files - ##ss_lst <- lapply(csvfiles, function(csv) read.csv(csv, header = TRUE, skip = 10, comment.char = '#')) - # use the first CSV file name as model name - m_name <- sub("(_\\d+)*$", '', sub("\\.[^.]*$", "", basename(csvfiles[1]))) - - sdate <- do.call(max, lapply(csvfiles, function(csv) file.info(csv)$mtime)) - sdate <- format(sdate, "%a %b %d %X %Y") # same format as date() - - chains <- length(ss_lst) - fnames <- names(ss_lst[[1]]) - n_save <- nrow(ss_lst[[1]]) - paridx <- paridx_fun(fnames) - lp__idx <- attr(paridx, 'meta')["lp__"] - par_fnames <- c(fnames[paridx], "lp__") - pars_oi <- unique_par(par_fnames) - dims_oi <- lapply(pars_oi, - function(i) { - pat <- paste('^', i, '(\\.\\d+)*$', sep = '') - i_fnames <- par_fnames[grepl(pat, par_fnames)] - get_dims_from_fnames(i_fnames, i) - }) - names(dims_oi) <- pars_oi - midx <- if (!col_major) multi_idx_row2colm(dims_oi) else 1:length(par_fnames) - if (chains > 1) { - if (!all(sapply(ss_lst[-1], function(i) identical(names(i), fnames)))) - stop('the CSV files do not have same parameters') - if (!all(sapply(ss_lst[-1], function(i) identical(length(i[[1]]), n_save)))) - stop('the number of iterations are not the same in all CSV files') - } - mode <- 0L - - samples <- lapply(ss_lst, - function(df) { - ss <- df[c(paridx, lp__idx)[midx]] - attr(ss, "sampler_params") <- df[setdiff(attr(paridx, 'meta'), lp__idx)] - ss - }) - par_fnames <- par_fnames[midx] - for (i in seq_along(samples)) { - attr(samples[[i]], "adaptation_info") <- cs_lst2[[i]]$adaptation_info - attr(samples[[i]], "args") <- - list(sampler_t = cs_lst2[[i]]$sampler_t, - chain_id = cs_lst2[[i]]$chain_id) - if (cs_lst2[[i]]$has_time) - attr(samples[[i]], "elapsed_time") <- get_time_from_csv(cs_lst2[[i]]$time_info) - } - - save_warmup <- sapply(cs_lst2, function(i) i$save_warmup) - warmup <- sapply(cs_lst2, function(i) i$warmup) - thin <- sapply(cs_lst2, function(i) i$thin) - iter <- sapply(cs_lst2, function(i) i$iter) - if (!all_int_eq(warmup) || !all_int_eq(thin) || !all_int_eq(iter)) - stop("not all iter/warmups/thin are the same in all CSV files") - n_kept0 <- 1 + (iter - warmup - 1) %/% thin - warmup2 <- 0 - if (max(save_warmup) == 0L) { # all equal to 0L - n_kept <- n_save - } else if (min(save_warmup) == 1L) { # all equals to 1L - warmup2 <- 1 + (warmup[1] - 1) %/% thin[1] - n_kept <- n_save - warmup2 - } - - if (n_kept0[1] != n_kept) { - warning("the number of iterations after warmup found (", n_kept, - ") does not match iter/warmup/thin from CSV comments (", - paste(n_kept0, collapse = ','), ")") - - if (n_kept < 0) { - warmup <- warmup + n_kept - n_kept <- 0 - mode <- 2L - } - n_kept0 <- n_save - iter <- n_save - - for (i in 1:length(cs_lst2)) { - cs_lst2[[i]]$warmup <- warmup - cs_lst2[[i]]$iter <- iter - } - } - - idx_kept <- if (warmup2 == 0) 1:n_kept else -(1:warmup2) - for (i in seq_along(samples)) { - m <- apply(samples[[i]][idx_kept,], 2, mean) - attr(samples[[i]], "mean_pars") <- m[-length(m)] - attr(samples[[i]], "mean_lp__") <- m["lp__"] - } - - perm_lst <- lapply(1:chains, function(id) sample.int(n_kept)) - - unique_par_fnames <- unique_par(par_fnames) - extract_draws <- function(warmup = TRUE) { - if(warmup) keep <- 1:warmup2 - else keep <- -c(1:warmup2) - draws <- sapply(unique_par_fnames, simplify = FALSE, FUN = function(u) { - mark <- grep(paste0("^", u, "\\."), colnames(samples[[1]]), value = TRUE) - out <- sapply(samples, FUN = function(s) { - if(u %in% colnames(s)) return(s[keep,u]) - return(as.matrix(s[keep,mark])) - }) - dims <- dim(out) - if(length(mark) == 0) { - out <- t(out) - dim(out) <- c(1L, rev(dims)) - dimnames(out) <- list(NULL, chains = NULL, iterations = NULL) - return(new("StanReal", name = u, theta = out, - type = "constrained parameter")) - } - dim(out) <- c(get_dims_from_fnames(mark, u), warmup2, length(samples)) - dims <- dim(out) - len <- length(dim(out)) - out <- aperm(out, perm = c(1:(len - 2L), len, len-1L)) - dimnames(out) <- vector("list", len) - names(dimnames(out))[len] <- "iterations" - names(dimnames(out))[len-1L] <- "chains" - return(new("StanMatrix", name = u, theta = out, - type = "constrained parameter")) - }) - draws <- c(draws, - sapply(colnames(attributes(samples[[1]])$sampler_params), - simplify = FALSE, FUN = function(u) { - out <- t(sapply(samples, function(s) { - return(attributes(s)$sampler_params[keep,u]) - })) - dim(out) <- c(1L, dim(out)) - dimnames(out) <- list(NULL, chains = NULL, - iterations = NULL) - return(new("StanReal", name = u, - theta = out, type = "diagnostic")) - })) - return(draws) - } - - sample_draws <- extract_draws(warmup = FALSE) - if(warmup2 > 0) warmup_draws <- extract_draws(warmup = TRUE) - else { - warmup_draws <- vector("list", length(sample_draws)) - names(warmup_draws) <- names(sample_draws) - } - - mass_matrix <- sapply(samples, FUN = function(s) { - m <- attributes(s)$adaptation_info - m <- sub("^#.*# ", "", m) - scan(text = m, sep = ",", quiet = TRUE) - }) - - seed <- as.integer(sapply(cs_lst2, FUN = function(x) x$seed)) - - return(new("StanFitMCMC", warmup_draws = warmup_draws, - sample_draws = sample_draws, seed = seed, - mass_matrix = mass_matrix)) - - # old output - sim = list(samples = samples, - iter = iter[1], - thin = thin[1], - warmup = warmup[1], - chains = chains, - n_save = rep(n_save, chains), - warmup2 = rep(warmup2, chains), - permutation = perm_lst, - pars_oi = pars_oi, - dims_oi = dims_oi, - fnames_oi = dotfnames_to_sqrfnames(par_fnames), - n_flatnames = length(par_fnames)) - null_dso <- new("cxxdso", sig = list(character(0)), dso_saved = FALSE, dso_filename = character(0), - modulename = character(0), system = R.version$system, cxxflags = character(0), - .CXXDSOMISC = new.env(parent = emptyenv())) - null_sm <- new("stanmodel", model_name = m_name, model_code = character(0), - model_cpp = list(), dso = null_dso) - - nfit <- new("stanfit", - model_name = m_name, - model_pars = pars_oi, - par_dims = dims_oi, - mode = mode, - sim = sim, - inits = list(), - stan_args = cs_lst2, - stanmodel = null_sm, - date = sdate, # not the time of sampling - .MISC = new.env(parent = emptyenv())) - return(nfit) -} diff --git a/rstan3/R/zzz.R b/rstan3/R/zzz.R deleted file mode 100644 index ea50ef932..000000000 --- a/rstan3/R/zzz.R +++ /dev/null @@ -1,22 +0,0 @@ -# This file is part of RStan -# Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017 Trustees of Columbia University -# -# RStan is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License -# as published by the Free Software Foundation; either version 3 -# of the License, or (at your option) any later version. -# -# RStan is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -.onLoad <- function(libname, pkgname) { } - -.onAttach <- function(...) { - packageStartupMessage("This is rstan3, a preview of version 3.x of rstan") -} -