R/libbi.R

In rbi: Interface to 'LibBi'

#' @rdname libbi
#' @name libbi
#' @title LibBi Wrapper
#' @description
#' \code{libbi} allows to call \code{LibBi}.
#' Upon creating a new libbi object, the following arguments can be given.
#' Once the instance is created, \code{LibBi} can be run through the
#'   \code{\link[rbi]{sample}}, \code{\link[rbi]{filter}}, or
#'   \code{\link{optimise}}, or \code{\link{rewrite}} methods. Note that
#'   \code{\link{libbi}} objects can be plotted using \code{\link{plot}} if the
#'   \code{rbi.helpers} package is loaded.
#'
#' @param model either a character vector giving the path to a model file
#'   (typically ending in ".bi"), or a \code{bi_model} object
#' @param path_to_libbi path to \code{LibBi} binary; by default it tries to
#'   locate the \code{libbi} binary using the \code{which} Unix command, after
#'   having loaded "~/.bashrc" if present; if unsuccessful it tries
#'   "~/PathToBiBin/libbi"; if unsuccessful again it fails.
#' @param dims any named dimensions, as list of character vectors
#' @param use_cache logical; whether to use the cache (default: true)
#' @param ... options passed to \code{\link{run.libbi}}
#' @return a new \code{\link{libbi}} object
#' @examples
#' bi_object <- libbi(model = system.file(package = "rbi", "PZ.bi"))
#' @seealso \code{\link{sample}}, \code{\link{filter}}, \code{\link{optimise}},
#'   \code{\link{rewrite}}
#' @export
libbi <- function(model, path_to_libbi, dims, use_cache = TRUE, ...) {
  libbi_dims <- list()
  if (!missing(dims)) {
    for (dim_name in names(dims)) {
      libbi_dims[[dim_name]] <- factor(dims[[dim_name]])
    }
  }

  if (missing(path_to_libbi)) path_to_libbi <- character(0)

  if (missing(model)) model <- bi_model()

  new_obj <- structure(list(
    options = list(),
    path_to_libbi = path_to_libbi,
    model = NULL,
    model_file_name = character(0),
    dims = libbi_dims,
    time_dim = character(0),
    coord_dims = list(),
    thin = 1,
    output_every = NA_real_,
    debug = FALSE,
    command = character(0),
    output_file_name = character(0),
    log_file_name = character(0),
    user_log_file = FALSE,
    timestamp = list(),
    run_flag = FALSE,
    error_flag = FALSE,
    use_cache = use_cache,
    supplement = NULL,
    .gc_env = emptyenv(),
    .cache = new.env(parent = emptyenv())
  ), class = "libbi")

  dot_options <- list(...)
  if ("norun" %in% names(dot_options) && dot_options[["norun"]]) {
    new_obj$model <- model
    return(new_obj)
  } else {
    return(do.call(
      run.libbi, c(
        list(x = new_obj, model = model, client = character(0)), dot_options
      )
    ))
  }
}

#' @export
run <- function(x, ...) UseMethod("run")
#' @rdname run
#' @name run
#' @title Using the LibBi wrapper to launch LibBi
#' @description
#' The method \code{run} launches \code{LibBi} with a particular set of command
#'   line arguments. Normally, this function would not be run by the user,
#'   but instead one of the client functions \code{\link{sample}},
#'   \code{\link{filter}}, or \code{\link{optimise}}, or \code{\link{rewrite}},
#'   which pass any options on to \code{run}. Note that any options specified
#'   here are stored in the \code{\link{libbi}} object and do not have to be
#'   specified again if another command is run on the object.
#'
#' @param x a \code{\link{libbi}} object; if this is not given, an empty
#'   \code{\link{libbi}} object will be created
#' @param client client to pass to LibBi
#' @param proposal proposal distribution to use; either "model" (default:
#'   proposal distribution in the model) or "prior" (propose from the prior
#'   distribution)
#' @param model either a character vector giving the path to a model file
#'   (typically ending in ".bi"), or a \code{bi_model} object; by default, will
#'   use any model given in \code{x}
#' @param fix any variable to fix, as a named vector
#' @param config path to a configuration file, containing multiple arguments
#' @param log_file_name path to a file to text file to report the output of
#'   \code{LibBi}; if set to an empty vector (\code{character(0)}) or an empty
#'   string (""), which is the default, a temporary log file will be generated
#' @param init initialisation of the model, either supplied as a list of values
#'   and/or data frames, or a (netcdf) file name, or a \code{\link{libbi}}
#'   object which has been run (in which case the output of that run is used).
#'   If the object given as \code{x} has been run before, it will be used here
#'   with \code{init-np} set to the last iteration of the previous run, unless
#'   \code{init} is given explicitly.
#' @param input input of the model, either supplied as a list of values and/or
#'   data frames, or a (netcdf) file name, or a \code{\link{libbi}} object which
#'   has been run (in which case the output of that run is used as input)
#' @param obs observations of the model, either supplied as a list of values
#'   and/or data frames, or a (netcdf) file name, or a \code{\link{libbi}}
#'   object which has been run (in which case the output of that run is used as
#'   observations)
#' @param time_dim The time dimension in any R objects that have been passed
#'   (\code{init}, \code{input}) and \code{obs}); if NULL (default), will be
#'   guessed from the given observation
#' @param coord_dims The coord dimension(s) in any \code{obs} R objects that
#'   have been passed; if NULL (default), will be guessed from the given
#'   observation file given
#' @param thin any thinning of MCMC chains (1 means all will be kept, 2 skips
#'   every other sample etc.); note that \code{LibBi} itself will write all data
#'   to the disk. Only when the results are read in with \code{\link{bi_read}}
#'   will thinning be applied.
#' @param output_every real; if given, \code{noutputs} will be set so that there
#'   is output every \code{output_every} time steps; if set to 0, only generate
#'   an output at the final time
#' @param chain logical; if set to TRUE and \code{x} has been run before, the
#'   previous output file will be used as \code{init} file, and \code{init-np}
#'   will be set to the last iteration of the previous run (unless
#'   target=="prediction"). This is useful for running inference chains.
#' @param seed Either a number (the seed to supply to \code{LibBi}), or a
#'   logical variable: TRUE if a seed is to be generated for \code{RBi}, FALSE
#'   if \code{LibBi} is to generate its own seed
#' @param debug logical; if TRUE, print more verbose messages and write all
#'   variables to the output file, irrespective of their setting of 'has_output'
#' @param ... list of additional arguments to pass to the call to \code{LibBi}.
#'   Any arguments starting with `enable`/`disable` can be specified as boolean
#'   (e.g., `assert=TRUE` or `cuda=TRUE`). Any `dry-` options can be specified
#'   with a `"dry"` argument, e.g., `dry="parse"`. Any options that would be
#'   specified with `with`/`without` can be specified as character vector to an
#'   option named `with`/`without`, respectively, e.g.
#'   with="transform-obs-to-state".
#' @seealso \code{\link{libbi}}
#' @return an updated \code{\link{libbi}} object, except if \code{client} is
#'   'rewrite', in which case invisible NULL will be returned but the rewritten
#'   model code printed
#' @importFrom ncdf4 nc_open nc_close ncvar_rename
#' @importFrom stats runif
#' @importFrom processx run
#' @export
run.libbi <- function(x, client, proposal = c("model", "prior"), model, fix,
                      config, log_file_name = character(0), init, input, obs,
                      time_dim = character(0), coord_dims = list(), thin,
                      output_every, chain = TRUE, seed = TRUE, debug = FALSE,
                      ...) {
  ## client options
  libbi_client_args <- list(
    sample = c(
      "target", "sampler", "nsamples", "nmoves", "tmoves",
      "sampler-resampler", "sample-ess-rel", "sample-stopper",
      "sample-stopper-threshold", "sample-stopper-max",
      "adapter", "adapter-scale", "adapter-ess-rel"
    ),
    optimise = c(
      "target", "optimiser", "simplex-size-real",
      "stop-size", "stop-steps"
    ),
    filter = c(
      "start-time", "end-time", "noutputs",
      "with-output-at-obs", "filter", "nparticles", "ess-rel",
      "resampler", "nbridges", "stopper", "stopper-threshold",
      "stopper-max", "stopper-block"
    ),
    rewrite = c()
  )

  all_client_args <- unique(unname(unlist(libbi_client_args)))

  ## both sample and optimise inherit from filter
  libbi_client_args[["sample"]] <- unique(c(
    libbi_client_args[["sample"]], libbi_client_args[["filter"]]
  ))
  libbi_client_args[["optimise"]] <- unique(c(
    libbi_client_args[["optimise"]], libbi_client_args[["filter"]]
  ))

  if (!missing(x) && "bi_model" %in% class(x)) {
    if (missing(model)) {
      model <- x
    } else {
      stop("Two 'bi_model' objects given as 'x' and 'model'.")
    }
  }
  if (missing(x) || ("bi_model" %in% class(x))) {
    x <- libbi(norun = TRUE)
  }

  ## assign stored arguments
  if (!missing(thin)) {
    x$thin <- thin
  }
  if (!missing(output_every)) {
    x$output_every <- output_every
  }
  if (!missing(debug)) {
    x$debug <- debug
  }

  new_options <- list(...)

  if (missing(config) || config == "") {
    config_file_options <- list()
  } else {
    config_str <- config
    if (substr(config_str, 1, 1) == "@") {
      config_str <- substr(config_str, 2, nchar(config_str))
    }
    if (file.exists(config_str)) {
      config_file_options <- paste(readLines(config_str), collapse = " ")
    } else {
      stop("Could not find config file ", config_str)
    }
  }

  libbi_seed <- integer(0)
  if (is.logical(seed)) {
    if (seed == TRUE) {
      libbi_seed <- ceiling(runif(1, -1, .Machine$integer.max - 1))
    }
  } else {
    libbi_seed <- seed
  }
  if (length(libbi_seed) > 0) new_options[["seed"]] <- libbi_seed

  if (is.null(getOption("libbi_args"))) {
    libbi_option_args <- list()
  } else {
    libbi_option_args <- getOption("libbi_args")
  }

  all_options <- option_list(
    libbi_option_args, config_file_options, x$options, new_options
  )

  ## check if 'model-file' is contained in any options
  if ("model-file" %in% names(all_options)) {
    if (is_empty(model)) {
      model <- bi_model(all_options[["model-file"]])
      warning(
        "A model has been passed via 'model-file'. It will be stored in the ",
        "libbi object internally. If the file passed as 'model-file' is ",
        "changed, 'model-file' needs to be passed again."
      )
    } else {
      warning(
        "'model-file' and 'model' options both provided. Will ignore ",
        "'model-file'."
      )
    }
  }

  ## get model
  if (!missing(model)) x$model <- model

  if (!("bi_model" %in% class(x$model))) {
    x$model <- bi_model(filename = x$model)
  }

  if (is_empty(x$model)) {
    stop("No model specified")
  }

  if ("build-dir" %in% names(all_options)) {
    x$options[["build-dir"]] <- all_options[["build-dir"]]
  } else {
    x <- create_working_folder(x)
  }

  x$model_file_name <- file.path(
    absolute_path(x$options[["build-dir"]]),
    paste0(get_name(x$model), ".bi")
  )

  args <- match.call()

  if (!is.na(x$output_every)) {
    if ("noutputs" %in% names(args)) {
      if (missing(output_every)) {
        x$output_every <- NA_real_
      } else {
        stop(
          "noutputs' and 'output_every' both given. Must be one or the other"
        )
      }
    } else {
      start_time <- ifelse(
        "start-time" %in% names(all_options), all_options[["start-time"]], 0
      )
      end_time <- ifelse(
        "end-time" %in% names(all_options), all_options[["end-time"]], 0
      )
      new_options[["noutputs"]] <- ifelse(
        x$output_every == 0, 0, (end_time - start_time) / x$output_every
      )
    }
  }

  ## read file options: input, init, obs
  file_types <- c("input", "init", "obs")
  file_args <- intersect(names(args), file_types)
  ## assign file args to global options
  for (arg in file_args) x$options[[arg]] <- get(arg)

  if (x$run_flag && length(x$output_file_name) == 1 &&
    file.exists(x$output_file_name)) {
    added_options <- option_list(new_options)
    init_file_given <-
      ("init" %in% file_args && !is.null(x$options[["init"]])) ||
        "init-file" %in% names(added_options)
    init_np_given <- "init-np" %in% names(added_options)
    init_given <- init_file_given || init_np_given
    if (missing(chain)) { ## if chain not specified, only chain if no init
      ## option is given
      chain <- !init_given
    }
    if (chain) {
      if (init_file_given) {
        warning(
          "init file given and 'chain = TRUE'. Will ignore 'init' option. ",
          "To use the 'init' option, set 'chain = FALSE'."
        )
      }
      if (init_np_given) {
        warning(
          "'init-np' given as new option and 'chain = TRUE'. Will ignore ",
          "'init-np' option. To use the 'init-np' option, set 'chain = FALSE'"
        )
      }
      if ("target" %in% names(all_options) &&
        all_options[["target"]] == "prediction") {
        read_init <- bi_read(x, type = c("param", "state", "obs"))
        np_dims <- bi_dim_len(x$output_file_name, "np")
        x$options[["nsamples"]] <- floor(np_dims / x$thin)
        if (x$thin > 1) {
          x$thin <- 1
        }
        x$options[["init"]] <- read_init
      } else {
        types <- "param"
        chain_init <-
          ("with-transform-initial-to-param" %in% names(all_options))
        if (chain_init) {
          types <- c(types, "state")
        }
        read_init <- bi_read(x, type = types, init_to_param = chain_init)
        ## only take last sample
        x$options[["init"]] <- extract_sample(read_init, "last")
      }
      file_args <- union(file_args, "init")
    }
  }

  if (!missing(time_dim)) {
    x$time_dim <- time_dim
  }
  if (!missing(coord_dims)) {
    x$coord_dims <- coord_dims
  }

  ## loop over global options that are file args
  for (file in file_args) {
    arg <- x$options[[file]]
    ## unset global option (we set the file option instead later)
    x$options[[file]] <- NULL
    ## attach data here
    x <- attach_data(x, file, arg, time_dim = time_dim, coord_dims = coord_dims)
  }

  all_options <- option_list(
    libbi_option_args, config_file_options,
    x$options, new_options
  )

  if (length(log_file_name) == 1 && log_file_name == "") {
    log_file_name <- character(0)
  }

  if (!x$user_log_file && length(log_file_name) == 0) {
    x$log_file_name <- tempfile(
      pattern = "output", fileext = ".txt",
      tmpdir = absolute_path(all_options[["build-dir"]])
    )
    x$user_log_file <- FALSE
  } else if (length(log_file_name) > 0) {
    x$log_file_name <- absolute_path(
      filename = log_file_name, dirname = getwd()
    )
    x$user_log_file <- TRUE
  }

  if (length(client) > 0) {
    ## re-read options
    ## remove arguments of other clients
    retain_options <- setdiff(
      names(all_options),
      setdiff(all_client_args, libbi_client_args[[client]])
    )
    all_options <- all_options[retain_options]

    ## adjust options
    if (client != "rewrite") {
      ## clear cache
      x$.cache <- new.env(parent = emptyenv())

      if (!("output-file" %in% names(all_options))) {
        x$output_file_name <- tempfile(
          pattern = paste(get_name(x$model), "output", sep = "_"),
          fileext = ".nc", tmpdir = absolute_path(all_options[["build-dir"]])
        )
      } else {
        x$output_file_name <- absolute_path(
          all_options[["output-file"]], getwd()
        )
      }
    }

    ## save options
    x$options <- all_options
    if (length(x$output_file_name) > 0 && nchar(x$output_file_name) > 0) {
      all_options[["output-file"]] <- x$output_file_name
    } else {
      all_options[["output-file"]] <- NULL
    }

    save_model <- x$model

    if (x$debug) x$model <- enable_outputs(x$model)

    proposal <- match.arg(proposal)
    if (proposal == "prior") {
      x$model <- propose_prior(x$model)
    }

    if (!missing(fix)) {
      x$model <- do.call(fix.bi_model, c(list(x = x$model), as.list(fix)))
    }

    write_model(x)
    all_options[["model-file"]] <- x$model_file_name

    if (x$debug) all_options[["verbose"]] <- TRUE
    verbose <-
      ("verbose" %in% names(all_options) && all_options[["verbose"]] == TRUE)

    run_args <- option_string(all_options)

    x$path_to_libbi <- locate_libbi(x$path_to_libbi)
    con <- file(
      ifelse(length(x$log_file_name) == 0, "", x$log_file_name), open = "w+"
    )
    cb_stdout <- function(line, proc) {
      if (x$debug || client == "rewrite") {
        cat(line, "\n")
      }
      writeLines(line, con)
      flush(con)
    }
    cb_stderr <- function(line, proc) {
      if (x$debug || (verbose && grepl("(^[0-9]+:|\\.\\.\\.$)", line))) {
        cat(line, "\n")
      }
      writeLines(line, con)
      flush(con)
    }
    x$command <- paste(x$path_to_libbi, client, paste(run_args, collapse = " "))
    if (verbose) {
      cat("Launching LibBi...", "\n")
    }
    if (!(client == "rewrite" && !debug)) {
      cb_stdout(x$command)
    }
    p <- tryCatch(
      processx::run(
        command = x$path_to_libbi, args = c(client, run_args),
        error_on_status = FALSE, wd = all_options[["build-dir"]],
        spinner = verbose,
        stdout_line_callback = cb_stdout,
        stderr_line_callback = cb_stderr
      ),
      finally = close(con)
    )
    if (p$status != 0) {
      error_lines <- strsplit(p$stderr, "\n")[[1]]
      error_msg <-
        sub(
          "^[[:space:]]+", "", ## remove leading spaces
          sub(
            "\\n$", "",
            grep(
              "(configure: error|Error| failed\\.$)", error_lines, value = TRUE
            )
          )
        )
      stop_msg <- paste0("LibBi terminated with \"", error_msg[1], "\".")
      if (length(x$log_file_name) > 0) {
        stop_msg <- paste0(
          stop_msg, "\nYou can view a full log using \"print_log('",
          x$log_file_name, "')\""
        )
      }
      stop(stop_msg)
    } else if (verbose) {
      cat("...LibBi has finished!", "\n")
    }
    x$error_flag <- FALSE

    ## don't store target
    x$options$target <- NULL

    if (client == "rewrite") {
      invisible(NULL)
    } else {
      x$run_flag <- TRUE
      x <- update(x)
      ## get original model back if it has been modified
      x$model <- save_model
      return(x)
    }
  } else {
    ## if run from the constructor, just write the model and add all the options
    write_model(x)
    x$options <- all_options
    return(x)
  }
}

#' @export
sample <- function(x, ...) UseMethod("sample")
#' @name sample
#' @rdname sample
#' @title Using the LibBi wrapper to sample
#' @description
#' The method \code{sample} launches \code{libbi} to sample from a (prior,
#'   posterior or joint) distribution. See the options to
#'   \code{\link{run.libbi}} for how to specify the various components of
#'   sampling with LibBi, and the LibBi manual for all options that can be
#'   passed when the client is \code{sample}.
#'
#' If \code{x} is given as a 'bi_model', a \code{\link{libbi}} object will be
#'   created from the model For the help page of the base R \code{sample}
#'   function, see \code{\link[base]{sample}}.
#' @param x a \code{\link{libbi}} or \code{\link{bi_model}} object, or the name
#'   of a file containing the model
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return an updated \code{\link{libbi}} object
#' @export
sample.libbi <- function(x, ...) {
  run.libbi(x, client = "sample", ...)
}
#' @rdname sample
#' @name sample
#' @export
sample.bi_model <- function(x, ...) {
  run.libbi(x, client = "sample", ...)
}
#' @export
sample.default <- function(x, ...) {
  base::sample(x, ...)
}

#' @export
filter <- function(x, ...) UseMethod("filter")
#' @rdname filter
#' @name filter
#' @title Using the LibBi wrapper to filter
#' @description
#' The method \code{filter} launches \code{libbi} to filter state trajectories.
#'   See the options to \code{\link{run.libbi}} for how to specify the various
#'   components of sampling with LibBi, and the LibBi manual for all options
#'   that can be passed when the client is \code{filter}.
#'
#' If \code{x} is given as a 'bi_model', a \code{\link{libbi}} object will be
#'   created from the model For the help page of the base R \code{filter}
#'   function, see \code{\link[stats]{filter}}.
#' @param x a \code{\link{libbi}} or \code{\link{bi_model}} object, or the name
#'   of a file containing the model
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return an updated \code{\link{libbi}} object
#' @export
filter.libbi <- function(x, ...) {
  run.libbi(x, client = "filter", ...)
}
#' @rdname filter
#' @name filter
#' @export
filter.bi_model <- function(x, ...) {
  run.libbi(x, client = "filter", ...)
}
#' @export
filter.default <- function(x, ...) {
  stats::filter(x, ...)
}

#' @export
optimise <- function(x, ...) UseMethod("optimise")
#' @rdname optimise
#' @name optimise
#' @title Using the LibBi wrapper to optimise
#' @description
#' The method \code{optimise} launches \code{libbi} to optimise the parameters
#'   with respect to the likelihood or posterior distribution. See the options
#'   to \code{\link{run.libbi}} for how to specify the various components of
#'   sampling with LibBi, and the LibBi manual for all options that can be
#'   passed when the client is \code{optimise}.
#'
#' If \code{x} is given as a 'bi_model', a \code{\link{libbi}} object will be
#'   created from the model For the help page of the base R \code{optimise}
#'   function, see \code{\link[stats]{optimise}}.
#' @param x a \code{\link{libbi}} or \code{link{bi_model}} object, or the name
#'   of a file containing the model
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return an updated \code{\link{libbi}} object
#' @export
optimise.libbi <- function(x, ...) {
  run.libbi(x, client = "optimise", ...)
}
#' @rdname optimise
#' @name optimise
#' @export
optimise.bi_model <- function(x, ...) {
  run.libbi(x, client = "optimise", ...)
}
#' @export
optimise.default <- function(x, ...) {
  stats::optimise(x, ...)
}

#' @export
rewrite <- function(x, ...) UseMethod("rewrite")
#' @rdname rewrite
#' @name rewrite
#' @title Using the LibBi wrapper to rewrite
#' @description
#' The method \code{rewrite} launches \code{LibBi} to rewrite a model to inspect
#'   its internal representation in \code{LibBi}
#'
#' If \code{x} is given as a 'bi_model', a \code{\link{libbi}} object will be
#'   created from the model
#' @param x a \code{\link{libbi}} or \code{\link{bi_model}} object, or the name
#'   of a file containing the model
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return a re-written \code{\link{bi_model}} object
#' @export
rewrite.libbi <- function(x, ...) {
  run.libbi(x, client = "rewrite", ...)
}
#' @rdname rewrite
#' @name rewrite
#' @export
rewrite.bi_model <- function(x, ...) {
  run.libbi(x, client = "rewrite", ...)
}

#' @export
simulate <- function(x, ...) UseMethod("simulate")
#' @rdname simulate
#' @name simulate
#' @title Using the LibBi wrapper to simulate
#' @description
#' The method \code{simulate} launches \code{LibBi} to simulate a model by
#'   passing `target="joint"` to \code{LibBi}
#'
#' If \code{x} is given as a 'bi_model', a \code{\link{libbi}} object will be
#'   created from the model
#' @param x a \code{\link{libbi}} or \code{\link{bi_model}} object, or the name
#'   of a file containing the model
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return an updated \code{\link{bi_model}} object
#' @export
simulate.libbi <- function(x, ...) {
  run.libbi(x, client = "sample", target = "joint", ...)
}
#' @rdname simulate
#' @name simulate
#' @export
simulate.bi_model <- function(x, ...) {
  run.libbi(x, client = "sample", target = "joint", ...)
}

#' @export
attach_data <- function(x, ...) UseMethod("attach_data")
#' @name attach_data
#' @rdname attach_data
#' @title Attach a new file or data set to a \code{\link{libbi}} object
#' @description
#' Adds an (output, obs, etc.) file to a \code{\link{libbi}} object. This is
#'   useful to recreate a \code{\link{libbi}} object from the model and output
#'   files of a previous run
#'
#' The \code{\link{bi_write}} options \code{append} and \code{overwrite}
#'   determine what exactly the file will contain at the end of this. If they
#'   are both FALSE (the default), any existing file will be ignored. If
#'   \code{append} is TRUE, the existing data in the file will be preserved, and
#'   any data set passed as \code{data} and not already in the file will be
#'   added. If \code{overwrite} is TRUE, existing data in the file will be
#'   preserved except for variables that exist in the passed \code{data}.
#' @param x a \code{\link{libbi}} object
#' @param file the type of the file to attach, one of "output", "obs", "input"
#'   or "init"
#' @param data name of the file to attach, or a list of data frames
#'   that contain the outputs; it will be assumed that this is already thinned
#' @param in_place if TRUE, replace the file in place if it already exists in
#'   the libbi object; this can speed up the operation if append=TRUE as
#'   otherwise the file will have to be read and used again; it should be used
#'   with care, though, as it can render existing \code{\link{libbi}} objects
#'   invalid as the files they are pointing to are changed.
#' @param quiet if TRUE, will suppress the warning message normally given if
#'   replace=TRUE and the file exists already
#' @param ... any options to \code{\link{bi_write}} (e.g., 'time_dim')
#' @inheritParams bi_write
#' @return an updated \code{\link{libbi}} object
#' @examples
#' bi <- libbi(model = system.file(package = "rbi", "PZ.bi"))
#' example_output <- bi_read(system.file(package = "rbi", "example_output.nc"))
#' bi <- attach_data(bi, "output", example_output)
#' @export
attach_data.libbi <- function(x, file, data, in_place = FALSE, append = FALSE,
                              overwrite = FALSE, quiet = FALSE,
                              time_dim = character(0),
                              coord_dims = list(), ...) {
  if (file == "output") {
    existing_file_name <- x[["output_file_name"]]
  } else {
    existing_file_name <- x[["options"]][[paste0(file, "-file")]]
  }

  if (in_place && length(existing_file_name) == 1) {
    if (!quiet) {
      warning(
        "attaching data with 'in_place=TRUE'. This should be used carefully, ",
        "as it may render an existing object pointing to the file invalid. ",
        "To suppress this message, use quiet = TRUE."
      )
    }
    target_file_name <- existing_file_name
  } else {
    target_file_name <-
      tempfile(
        pattern = paste(get_name(x$model), file, sep = "_"),
        fileext = ".nc", tmpdir = absolute_path(x$options[["build-dir"]])
      )
    if (length(existing_file_name) == 1 && (append || overwrite)) {
      file.copy(existing_file_name, target_file_name)
    }
  }

  if (is.character(data) || ("libbi" %in% class(data))) {
    if (is.character(data)) {
      source_file_name <- data
    } else {
      tryCatch(assert_files(data),
        error = function(e) stop("Error adding ", file, " file\n", e)
      )
      if (length(data$time_dim) == 0 && length(time_dim) == 0) {
        time_dim <- x$time_dim
      } else {
        if (length(time_dim) == 0) {
          time_dim <- data$time_dim
        }
      }

      if (length(time_dim) > 0) {
        if (length(x$time_dim) > 0 && x$time_dim != time_dim) {
          warning(
            "Given time dimension will override the time dimension ",
            "of the object it is being added to."
          )
        }
        x$time_dim <- data$time_dim
      }

      if (file %in% c("obs", "input")) {
        if (length(coord_dims) == 0 && length(data$coord_dims) == 0) {
          coord_dims <- x$coord_dims
        } else {
          if (length(coord_dims) == 0) {
            coord_dims <- data$coord_dims
          }
        }
      }

      if (length(coord_dims) > 0) {
        for (coord_dim in names(coord_dims)) {
          if (!is.null(x$coord_dims[[coord_dim]]) &&
            any(x$coord_dims[[coord_dim]] != coord_dims[[coord_dim]])) {
            warning(
              "Given coord dimension ", coord_dim,
              " will override a coord dimension of the same name in",
              " passed libbi object"
            )
          }
        }
        x$coord_dims <- coord_dims
      }
      source_file_name <- data$output_file_name
    }

    if (file %in% c("obs", "input")) {
      vars <- bi_read(
        data, vars = var_names(x$model, type = "obs"), coord_dims = coord_dims
      )
    } else if (append || overwrite) {
      vars <- bi_read(data)
    } else {
      file.copy(source_file_name, target_file_name)
    }
  } else if ("list" %in% class(data) && length(data) > 0) {
    vars <- data
  } else { ## can't do anything with given 'data'
    if (file == "output") {
      x$run_flag <- FALSE
      x$output_file_name <- character(0)
      x$timestamp[["output"]] <- NULL
      x$.cache <- new.env(parent = emptyenv())
    } else {
      x$timestamp[[file]] <- NULL
      x$options[[paste(file, "file", sep = "-")]] <- NULL
    }
    return(x)
  }

  if ((append || overwrite || "list" %in% class(data) ||
       file %in% c("obs", "input")) &&
    length(vars) > 0) {
    write_opts <- list(filename = target_file_name, variables = vars)
    if (length(x$time_dim) == 0) {
      write_opts[["guess_time"]] <- TRUE
    } else {
      write_opts[["time_dim"]] <- x$time_dim
    }
    if (file %in% c("obs", "input")) {
      if (length(x$coord_dim) > 0) {
        write_opts[["coord_dim"]] <- x$coord_dim
      }
    }
    write_opts[["dim_factors"]] <- x$dims
    write_opts[["append"]] <- append
    write_opts[["overwrite"]] <- overwrite
    ## add any options passed
    added_options <- list(...)
    for (option in names(added_options)) {
      write_opts[[option]] <- added_options[[option]]
    }
    file_dims <- do.call(bi_write, write_opts)
    x$dims[names(file_dims$dims)] <- file_dims$dims
    if (!is.null(file_dims$time_dim)) x$time_dim <- file_dims$time_dim
    if (!is.null(file_dims$coord_dims)) {
      x$coord_dims[names(file_dims$coord_dims)] <- file_dims$coord_dims
    }
  }

  if (file == "output") {
    x$output_file_name <- target_file_name
    x$run_flag <- TRUE
    x$timestamp[["output"]] <- file.mtime(x$output_file_name)
    x$.cache <- new.env(parent = emptyenv())
    x$thin <- 1 ## output file assumed to be already be thinned
  } else {
    x$options[[paste0(file, "-file")]] <- target_file_name
    x$timestamp[[file]] <- file.mtime(target_file_name)
  }
  return(x)
}

#' @export
save_libbi <- function(x, ...) UseMethod("save_libbi")
#' @name save_libbi
#' @rdname save_libbi
#' @title Write results of a \code{LibBi} run to an RDS file
#' @description
#' This saves all options, files and outputs of a \code{LibBi} run to an RDS
#'   file specified
#'
#' @param x a \code{\link{libbi}} object
#' @param name name of the RDS file(s) to save to. If \code{split=TRUE}, this
#'   will be taken as a base for the names of the files to be created, e.g.
#'   'dir/name.rds' to create files of the form name_....rds in directory 'dir'.
#' @param supplement any supplementary data to save
#' @param split Logical, defaults to \code{FALSE}. Should the objects from the
#'  \code{LibBi} run be saved separately in a folder.
#' @param ... any options to \code{\link{saveRDS}}
#' @return the return value of \code{saveRDS}, i.e. NULL invisibly
#' @export
save_libbi.libbi <- function(x, name, supplement, split = FALSE, ...) {
  if (missing(name)) {
    stop("Need to specify a name")
  }

  assert_files(x)

  if (split) {
    folder <- dirname(name)
    file_base <- sub("\\.rds$", "", basename(name))

    if (!dir.exists(folder)) {
      stop(
        "The folder specified for saving the Libbi object into does not exist."
      )
    }
  }

  save_obj <- list(
    model = x$model,
    internals = list(
      dims = x$dims,
      time_dim = x$time_dim,
      coord_dims = x$coord_dims,
      output_every = x$output_every,
      debug = x$debug
    ),
    supplement = x$supplement,
    output = bi_read(x)
  )

  options <- x$options

  for (file_type in c("init", "input", "obs")) {
    file_option <- paste(file_type, "file", sep = "-")
    if (file_option %in% names(x$options)) {
      save_obj[[file_type]] <- bi_read(x, file = file_type)
      options[[file_option]] <- NULL
    }
  }

  options[["build-dir"]] <- NULL
  save_obj[["options"]] <- options

  save_obj[["log"]] <- character(0)
  if (length(x$log_file_name) > 0) {
    if (file.exists(x$log_file_name)) {
      save_obj[["log"]] <- readLines(x$log_file_name)
    }
  }

  if (!missing(supplement)) save_obj[["supplement"]] <- supplement

  if (split) {
    filenames <- list()
    for (i in names(save_obj)) {
      if (i == "output") {
        filenames[[i]] <- list()
        for (j in names(save_obj[[i]])) {
          filenames[[i]][[j]] <- paste0(
            paste(file_base, i, j, sep = "_"), ".rds"
          )
          saveRDS(
            save_obj[[i]][[j]], file.path(folder, filenames[[i]][[j]]), ...
          )
        }
      } else {
        filenames[[i]] <- paste0(paste(file_base, i, sep = "_"), ".rds")
        saveRDS(save_obj[[i]], file.path(folder, filenames[[i]]), ...)
      }
    }
    saveRDS(list(filenames = filenames, split = c()), name, ...)
  } else {
    saveRDS(save_obj, name, ...)
  }
}

#' @rdname read_libbi
#' @name read_libbi
#' @title Read results of a \code{LibBi} run from an RDS file or from a folder.
#'   This completely reconstructs the saved \code{LibBi} object
#' @description
#' This reads all options, files and outputs of a \code{LibBi} run from a
#'   specified RDS file or folder (if \code{split = TRUE} has been used with
#'   \code{save_libbi}).
#'
#' @param name name of the RDS file(s) to read
#' @param ... any extra options to pass to \code{\link{libbi}} when creating the
#'   new object
#' @return a new \code{\link{libbi}} object
#' @export
read_libbi <- function(name, ...) {
  if (missing(name)) {
    stop("Need to specify a file or folder to read from")
  }

  read_obj <- readRDS(name)

  if ("split" %in% names(read_obj)) { ## object has been split
    folder <- dirname(name)

    filenames <- read_obj[["filenames"]]

    read_obj <- list()
    for (obj_name in names(filenames)) {
      if (obj_name == "output") {
        for (var_name in names(filenames[["output"]])) {
          obj <- list(readRDS(file.path(
            folder, filenames[["output"]][[var_name]]
          )))
          names(obj) <- var_name
          read_obj[["output"]] <- c(read_obj[["output"]], obj)
        }
      } else {
        obj <- readRDS(file.path(folder, filenames[[obj_name]]))
        if (obj_name == "internals") {
          read_obj[names(obj)] <- obj
        } else {
          read_obj[[obj_name]] <- obj
        }
      }
    }
  }

  libbi_options <- list(...)

  pass_options <- c(
    "model", "dims", "time_dim", "coord_dims",
    "thin", "output-every", "init", "input", "obs"
  )

  for (option in pass_options) {
    if (!(option %in% names(libbi_options)) &&
      option %in% names(read_obj)) {
      libbi_options[[option]] <- read_obj[[option]]
    }
  }

  ## pass options directly
  if ("options" %in% names(read_obj)) {
    for (option in names(read_obj[["options"]])) {
      if (!(option %in% names(libbi_options))) {
        libbi_options[[option]] <- read_obj[["options"]][[option]]
      }
    }
  }

  libbi_options[["build-dir"]] <- NULL
  new_obj <- do.call(libbi, libbi_options)

  ## write output file
  new_obj <- attach_data(
    new_obj, file = "output", data = read_obj$output,
    time_dim = libbi_options$time_dim
  )

  ## write log file
  if ("log" %in% names(read_obj)) {
    writeLines(read_obj[["log"]], new_obj$log_file_name)
  }

  new_obj$supplement <- read_obj$supplement

  return(new_obj)
}

#' @export
#' @name print
#' @title Print information about a \code{\link{libbi}} object
#' @description
#' This prints the model name, basic information such as number of iterations
#'   and timesteps run, as well as a list of variables.
#' @param x a \code{\link{libbi}} object
#' @param verbose logical; if TRUE, locations of files and working folder should
#'   be printed
#' @param ... ignored
#' @rdname print
#' @return nothing (invisible NULL)
print.libbi <- function(x, verbose = FALSE, ...) {
  cat("Wrapper around LibBi\n")
  if (verbose) {
    cat("* path to working folder:", x$options[["build-dir"]], "\n")
    cat("* path to model file:", x$model_file_name, "\n")
    if (length(x$output_file_name) > 0) {
      cat("* path to output_file:", x$output_file_name, "\n")
    }
  }
  cat("======================\n")
  cat("Model: ", get_name(x$model), "\n")
  if (x$run_flag) {
    if (length(x$output_file_name) > 0 && file.exists(x$output_file_name)) {
      contents <- bi_contents(x$output_file_name)
      if ("clock" %in% contents) {
        clock <- bi_read(x, "clock")[["clock"]]
        cat("Run time: ", clock / 1e6, " seconds\n")
      }
      states <- intersect(contents, var_names(x$model, type = "state"))
      noises <- intersect(contents, var_names(x$model, type = "noise"))
      params <- intersect(contents, var_names(x$model, type = "param"))
      obs <- intersect(contents, var_names(x$model, type = "obs"))
      niterations <- bi_dim_len(x$output_file_name, "np")
      if (niterations > 0) {
        cat("Number of samples: ", niterations, "\n")
      }
      if (length(states) > 0) {
        cat("State trajectories recorded: ", paste(states, sep = ", "), "\n")
      }
      if (length(noises) > 0) {
        cat("Noise trajectories recorded: ", paste(noises, sep = ", "), "\n")
      }
      if (length(obs) > 0) {
        cat("Observation trajectories recorded: ", paste(obs, sep = ", "), "\n")
      }
      if (length(params) > 0) {
        cat("Parameters recorded: ", paste(params), "\n")
      }
    } else {
      cat("* No output file\n")
    }
  } else {
    cat("* LibBi has not been run yet\n")
  }
}

#' @export
#' @name print_log
#' @title Print the log file a \code{\link{libbi}} object
#' @description
#' This is useful for diagnosis after a \code{\link{libbi}} run
#' @param x a \code{\link{libbi}} object, or the name of the log file of a
#'   \code{\link{libbi}} run.
#' @rdname print_log
#' @return nothing (invisible NULL)
print_log <- function(x) {
  if ("libbi" %in% class(x)) {
    file_name <- x$log_file_name
    if (!file.exists(file_name)) {
      stop("Log file '", x$log_file_name, " does not exist.")
    }
  } else if (is.character(x)) {
    if (file.exists(x)) {
      file_name <- x
    } else {
      stop(
        "If 'x' is a character object, it must point to an existing file. ",
        x, "does not exist."
      )
    }
  } else {
    stop("'x' must be a 'libbi' or 'character' object")
  }

  lines <- readLines(file_name)
  for (i in seq_along(lines)) {
    cat(lines[i], "\n")
  }
}

#' @name summary
#' @rdname summary
#' @title Print summary information about a \code{\link{libbi}} object
#' @description
#' This reads in the output file of the \code{\link{libbi}} object (which has
#'   been run before) and prints summary information of parameters
#' @param object a \code{\link{libbi}} object
#' @param type one of "param" (default), "state", "noise" or "obs", the variable
#'   type to summarise
#' @param quantiles quantiles to calculate (default: quartiles); minimum,
#'   median, mean and maximum are always calculated
#' @param na.rm logical; if true, any \code{na} and \code{NaN}'s are removed
#'   before calculations are performed
#' @param ... ignored
#' @importFrom data.table data.table setDF rbindlist
#' @importFrom stats median quantile
#' @export
#' @return nothing (invisible NULL)
summary.libbi <- function(object, type = c("param", "state", "noise", "obs"),
                          quantiles = c(0.25, 0.75), na.rm = FALSE, ...) { # nolint
  type <- match.arg(type)
  vars <- c(bi_read(object, type = type))

  value <- NULL ## to prevent data.table errors

  summary_table <- rbindlist(lapply(names(vars), function(var) {
    object <- vars[[var]]
    summarise_columns <- setdiff(colnames(object), c("np", "value"))
    dt <- data.table(object)[, c(
      list(
        var = var,
        `Min.` = min(value, na.rm = na.rm),
        `Mean` = mean(value, na.rm = na.rm),
        `Max.` = max(value, na.rm = na.rm)
      ),
      as.list(quantile(
        value,
        union(0.5, quantiles),
        na.rm = na.rm
      ))
    ), by = summarise_columns]
    return(dt)
  }))
  ## reorder table
  numeric_columns <- c(
    "Min.",
    paste0(quantiles[quantiles < 0.5] * 100, "%"),
    "50%", "Mean",
    paste0(quantiles[quantiles > 0.5] * 100, "%"),
    "Max."
  )
  summary_table <- summary_table[, c(
    "var",
    setdiff(
      colnames(summary_table),
      c("var", numeric_columns)
    ),
    numeric_columns
  ), with = FALSE]
  setnames(
    summary_table,
    c("25%", "50%", "75%"),
    c("1st Qu.", "Median", "3rd Qu."),
    skip_absent = TRUE
  )
  setDF(summary_table)

  return(summary_table)
}

#' @keywords internal
assert_files <- function(x, ...) UseMethod("assert_files")
#' @name assert_files
#' @rdname assert_files
#' @title Check that a LibBi wrapper has valid output
#' @description
#' This checks that the \code{\link{libbi}} object given has been run (via
#'   \code{\link{sample}}, \code{\link{filter}} or \code{\link{optimize}})) and
#'   the output file has not been modified since.
#' @param x a \code{\link{libbi}} object
#' @param ... ignored
#' @keywords internal
#' @return nothing (invisible NULL); an error will be thrown if there is a
#'   problem
assert_files.libbi <- function(x, ...) {
  if (!x$run_flag) {
    stop(
      "The libbi object must be run first (using sample, filter or optimise)."
    )
  }
  if (length(x$output_file_name) == 0 || !file.exists(x$output_file_name)) {
    stop("The libbi object does not contain an output file.")
  } else {
    if ("output" %in% names(x$timestamp) &&
      x$timestamp[["output"]] < file.mtime(x$output_file_name)) {
      stop(
        "Output file ", x$output_file_name,
        " has been modified since LibBi was run."
      )
    }
  }

  for (file_option in grep("-file$", names(x$options), value = TRUE)) {
    file_type <- sub("-file$", "", file_option)
    if (file.exists(x$options[[file_option]])) {
      if (file_type %in% names(x$timestamp) &&
        x$timestamp[[file_type]] < file.mtime(x$options[[file_option]])) {
        stop(
          file_type, " file ", x$options[[file_option]],
          " has been modified since LibBi was run. You can use",
          " the 'update' function to accept this change."
        )
      }
    } else {
      stop(file_type, " file '", x$options[[file_option]], " does not exist.")
    }
  }
}

#' @export
predict <- function(x, ...) UseMethod("predict")
#' @name predict
#' @rdname predict
#' @title Using the LibBi wrapper to predict
#' @description
#' The method \code{predict} is an alias for \code{sample(target="prediction")}.
#'   Usually, an \code{init} object or file should be given containing posterior
#'   samples.
#'
#' For the help page of the base R \code{optimise} function, see
#'   \code{\link[stats]{optimise}}.
#' @export
#' @param x a \code{\link{libbi}} object
#' @param ... any arguments to be passed to \code{\link{sample}}
#' @return an updated \code{\link{libbi}} object
predict.libbi <- function(x, ...) {
  sample(x, target = "prediction", ...)
}
#' @export
predict.default <- function(x, ...) {
  stats::predict(x, ...)
}


##' Sample observations from a LibBi model that has been run
##'
##' @param x a \code{\link{libbi}} object
##' @param ... any options to pass to LibBi
##' @return the original \code{\link{libbi}} object with added variables in the
##'   output file for sampled observations
##' @author Sebastian Funk
##' @export
sample_obs <- function(x, ...) {
  if (!("libbi" %in% class(x))) {
    stop("'x' must be a 'libbi' object")
  }
  out <- bi_read(x)
  out$clock <- NULL ## we don't want to overwrite the clock later
  ## remove transition
  sample_model <- remove_lines(x$model, "transition", preserve_shell = TRUE)
  ## turn outputs into inputs
  pr <- attach_data(x, file = "input", out, append = TRUE)
  ## predict
  pr <- predict(pr,
    model = sample_model,
    with = c("output-at-obs", "transform-obs-to-state"), ...
  )
  ## attach outputs back
  pr <- attach_data(pr, "output", data = out, append = TRUE)
  pr$options[["input-file"]] <- x$options[["input-file"]]
  pr$options[["with-output-at-obs"]] <- x$options[["with-output-at-obs"]]
  pr$options[["without-output-at-obs"]] <- x$options[["without-output-at-obs"]]
  pr$options[["with-transform-obs-to-state"]] <-
    x$options[["with-transform-obs-to-state"]]
  pr$options[["without-transform-obs-to-state"]] <-
    x$options[["without-transform-obs-to-state"]]
  pr$model <- x$model

  return(pr)
}

#' @export
join <- function(x, ...) UseMethod("join")
#' @name join
#' @rdname join
#' @title Join multiple \code{\link{libbi}} objects
#' @description
#' This function can be used to join multiple \code{\link{libbi}} objects into
#'   one (e.g., parallel MCMC runs into one long change)
#' @export
#' @param x a \code{\link{libbi}} object
#' @param ... ignored
#' @return an joined \code{\link{libbi}} object
join.libbi <- function(x, ...) {
  files_to_join <- list(...)
  output <- bi_read(x)
  for (to_join in files_to_join) {
    join_output <- bi_read(to_join)
    for (var in intersect(names(output), names(join_output))) {
      if (var == "clock") {
        output$clock <- output$clock + join_output$clock
      } else if (
        is.data.frame(output[[var]]) &&
          is.data.frame(join_output[[var]]) &&
          "np" %in% colnames(output[[var]]) &&
          "np" %in% colnames(join_output[[var]])
      ) {
        join_output[[var]]$np <- join_output[[var]]$np +
          max(output[[var]]$np) + 1
        output[[var]] <- rbind(output[[var]], join_output[[var]])
      }
    }
    for (var in setdiff(names(output), names(join_output))) {
      output[[var]] <- NULL
    }
  }
  attach_data(x, file = "output", output)
}

#' @rdname logLik
#' @name logLik
#' @title Using the LibBi wrapper to logLik
#' @description
#' The method \code{logLik} extracts the log-likelihood of a \code{libbi}
#'   object. This can be done, for example, after a call to
#'   \code{\link[rbi]{sample}} to inspect the chain log-likelihoods.
#'
#' For the help page of the base R \code{logLik} function, see
#'   \code{\link[stats]{logLik}}.
#' @param object a \code{\link{libbi}} object
#' @param ... options to be passed to \code{\link{run.libbi}}
#' @return a vector of log-likelihood
#' @export
logLik.libbi <- function(object, ...) {
  assert_files(object)
  res <- bi_read(object)
  if (is.vector(res$loglikelihood)) {
    return(res$loglikelihood)
  } else {
    return(res$loglikelihood$value)
  }
}

#' @export
update <- function(x, ...) UseMethod("update")
#' @rdname update
#' @name update
#' @title Update a libbi object
#' @description
#' This updates all the time stamps in a libbi object; it is useful after
#' (input, output, etc.) files have been changed outside the object itself.
#'
#' @param x a \code{\link{libbi}} object
#' @param ... ignored
#' @return a \code{\link{libbi}} object with updated timestamps
#' @export
update.libbi <- function(x, ...) {
  x$timestamp <- list()
  if ("options" %in% names(x)) {
    file_options <- setdiff(
      grep("-file$", names(x$options), value = TRUE), "output-file"
    )
    for (file_option in file_options) {
      if (file.exists(x$options[[file_option]])) {
        file_type <- sub("-file$", "", file_option)
        x$timestamp[[file_type]] <- file.mtime(x$options[[file_option]])
      }
    }
  }
  if (length(x$output_file_name) > 0) {
    if (file.exists(x$output_file_name)) {
      x$timestamp[["output"]] <- file.mtime(x$output_file_name)
    }
  }
  return(x)
}
#' @export
update.default <- function(x, ...) {
  stats::update(x, ...)
}


#' @title Internal function to create a temporary working folder
#' @description
#' The folder will be removed when the object is destroyed
#' @param x a \code{\link{libbi}} object
#' @return a \code{\link{libbi}} object with updated working folder
#' @keywords internal
create_working_folder <- function(x) {
  x$options[["build-dir"]] <- tempfile(pattern = paste(get_name(x$model)))
  dir.create(x$options[["build-dir"]])
  ## make sure temporary folder gets deleted upon garbage collection
  x$.gc_env <- new.env() ## dummy environment
  x$.gc_env$folder <- x$options[["build-dir"]]
  reg.finalizer(x$.gc_env, function(env) {
    unlink(env$folder, recursive = TRUE)
  }, onexit = TRUE)
  return(x)
}