R/C-generator.R

In SimInf: A Framework for Data-Driven Stochastic Disease Spread Simulations

## This file is part of SimInf, a framework for stochastic
## disease spread simulations.
##
## Copyright (C) 2015 Pavol Bauer
## Copyright (C) 2017 -- 2019 Robin Eriksson
## Copyright (C) 2015 -- 2019 Stefan Engblom
## Copyright (C) 2015 -- 2022 Stefan Widgren
##
## SimInf 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.
##
## SimInf 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, see <https://www.gnu.org/licenses/>.

##' Generate C code for a heading with timestamp and SimInf version
##' @return character vector with C code.
##' @importFrom utils packageVersion
##' @noRd
C_heading <- function() {
    c(sprintf("/* Generated by SimInf (v%s) */", packageVersion("SimInf")),
      "")
}

##' Generate C code with include directives
##' @return character vector with C code.
##' @noRd
C_include <- function() {
    c("#include <R_ext/Rdynload.h>",
      "#include \"SimInf.h\"",
      "")
}

##' Generate C code for definitions in the heading
##' @return character vector with C code.
##' @noRd
C_define <- function() {
    c("/**",
      " * Make sure the necessary macros are defined so that the",
      " * compiler can replace them when compiling the model.",
      " * 'SIMINF_MODEL_RUN' defines the function name of the function",
      " * that will be called from R to run a trajectory of the model.",
      " * 'SIMINF_R_INIT' is the name of the function that R will call",
      " * when this model is loaded into R. 'SIMINF_FORCE_SYMBOLS'",
      " * defines whether R allows the entry point for the run function",
      " * to be searched for as a character string.",
      " * If this file is compiled from SimInf (when calling run), the",
      " * macros are defined by SimInf before calling 'R CMD SHLIB'.",
      " * If this file is compiled as part of a package, then the",
      " * definitions are set in the variable 'PKG_CPPFLAGS' in",
      " * 'src/Makevars' and 'src/Makevars.in'.",
      " */",
      "#if !defined(SIMINF_MODEL_RUN)",
      "#  error Definition for 'SIMINF_MODEL_RUN' is missing.",
      "#endif",
      "#if !defined(SIMINF_R_INIT)",
      "#  error Definition for 'SIMINF_R_INIT' is missing.",
      "#endif",
      "#if !defined(SIMINF_FORCE_SYMBOLS)",
      "#  error Definition for 'SIMINF_FORCE_SYMBOLS' is missing.",
      "#endif",
      "")
}

##' Generate C code for the model transition rate functions
##'
##' @param transitions data for the transitions.
##' @return character vector with C code.
##' @noRd
C_trFun <- function(transitions) {
    parameters <- c("    const int *u,",
                    "    const double *v,",
                    "    const double *ldata,",
                    "    const double *gdata,",
                    "    double t)")

    lines <- character(0)
    for (i in seq_len(length(transitions))) {
        lines <- c(lines,
                   "/**",
                   " * @param u The compartment state vector in the node.",
                   " * @param v The continuous state vector in the node.",
                   " * @param ldata The local data vector in the node.",
                   " * @param gdata The global data vector.",
                   " * @param t Current time.",
                   " * @return propensity.",
                   " */",
                   sprintf("static double trFun%i(", i),
                   parameters,
                   "{",
                   sprintf("    return %s;", transitions[[i]]$propensity),
                   "}",
                   "")
    }
    lines
}

##' Generate C code for a SimInf model post-time-step function
##'
##' @param pts_fun optional character vector with C code for the post
##'     time step function. The C code should contain only the body of
##'     the function i.e. the code between the opening and closing
##'     curly brackets.
##' @return character vector with C code.
##' @noRd
C_ptsFun <- function(pts_fun) {
    if (is.null(pts_fun))
        pts_fun <- "    return 0;"

    if (!is.character(pts_fun))
        stop("'pts_fun' must be a character vector.", call. = FALSE)

    f <- textConnection(pts_fun)
    lines <- readLines(f)
    close(f)

    c("/**",
      " * Post time step function.",
      " *",
      " * @param v_new If a continuous state vector is used by a model,",
      " *        this is the new continuous state vector in the node after",
      " *        the post time step.",
      " * @param u The compartment state vector in the node.",
      " * @param v The current continuous state vector in the node.",
      " * @param ldata The local data vector in the node.",
      " * @param gdata The global data vector that is common to all nodes.",
      " * @param node The node index. Note the node index is zero-based,",
      " *        i.e., the first node is 0.",
      " * @param t Current time in the simulation.",
      " * @return error code (<0), or 1 if node needs to update the",
      " *         transition rates, or 0 when it doesn't need to update",
      " *         the transition rates.",
      " */",
      "static int ptsFun(",
      "    double *v_new,",
      "    const int *u,",
      "    const double *v,",
      "    const double *ldata,",
      "    const double *gdata,",
      "    int node,",
      "    double t)",
      "{",
           lines,
      "}",
      "")
}

##' Generate C code for a SimInf model run function
##'
##' @param transitions data for the transitions.
##' @return character vector with C code.
##' @noRd
C_run <- function(transitions) {
    c("/**",
      " * Run a trajectory of the model.",
      " *",
      " * @param model The model.",
      " * @param solver The name of the numerical solver.",
      " * @return A model with a trajectory attached to it.",
      " */",
      "static SEXP SIMINF_MODEL_RUN(SEXP model, SEXP solver)",
      "{",
      "    static SEXP(*SimInf_run)(SEXP, SEXP, TRFun*, PTSFun) = NULL;",
      sprintf("    TRFun tr_fun[] = {%s};",
              paste0("&trFun", seq_len(length(transitions)), collapse = ", ")),
      "",
      "    if (!SimInf_run) {",
      "        SimInf_run = (SEXP(*)(SEXP, SEXP, TRFun*, PTSFun))",
      "            R_GetCCallable(\"SimInf\", \"SimInf_run\");",
      "",
      "        if (!SimInf_run) {",
      "            Rf_error(\"Cannot find function 'SimInf_run'.\");",
      "        }",
      "    }",
      "",
      "    return SimInf_run(model, solver, tr_fun, &ptsFun);",
      "}",
      "")
}

##' Generate C code for the calldef for registering native routines
##' @return character vector with C code.
##' @noRd
C_calldef <- function() {
    c("/**",
      " * A NULL-terminated array of routines to register for the .Call",
      " * interface, see section '5.4 Registering native routines' in",
      " * the 'Writing R Extensions' manual.",
      " */",
      "static const R_CallMethodDef callMethods[] =",
      "{",
      "    SIMINF_CALLDEF(SIMINF_MODEL_RUN, 2),",
      "    {NULL, NULL, 0}",
      "};",
      "")
}

##' Generate C code for the R init function for registering native
##' routines
##'
##' @return character vector with C code.
##' @noRd
C_R_init <- function() {
    c("/**",
      " * This routine will be invoked when R loads the shared object/DLL,",
      " * see section '5.4 Registering native routines' in the",
      " * 'Writing R Extensions' manual.",
      " */",
      "void SIMINF_R_INIT(DllInfo *info)",
      "{",
      "    R_registerRoutines(info, NULL, callMethods, NULL, NULL);",
      "    R_useDynamicSymbols(info, FALSE);",
      "    R_forceSymbols(info, SIMINF_FORCE_SYMBOLS);",
      "}")
}

##' Generate C code for mparse
##'
##' @param transitions data for the transitions.
##' @param pts_fun optional character vector with C code for the post
##'     time step function. The C code should contain only the body of
##'     the function i.e. the code between the opening and closing
##'     curly brackets.
##' @return character vector with C code.
##' @noRd
C_code_mparse <- function(transitions, pts_fun) {
    c(C_heading(),
      C_include(),
      C_define(),
      C_trFun(transitions),
      C_ptsFun(pts_fun),
      C_run(transitions),
      C_calldef(),
      C_R_init())
}

##' Extract the C code from a \code{SimInf_model} object
##'
##' @param model The \code{SimInf_model} object to extract the C code
##'     from.
##' @return Character vector with C code for the model.
##' @export
##' @examples
##' ## Use the model parser to create a 'SimInf_model' object that
##' ## expresses an SIR model, where 'b' is the transmission rate and
##' ## 'g' is the recovery rate.
##' model <- mparse(transitions = c("S -> b*S*I/(S+I+R) -> I", "I -> g*I -> R"),
##'                 compartments = c("S", "I", "R"),
##'                 gdata = c(b = 0.16, g = 0.077),
##'                 u0 = data.frame(S = 99, I = 1, R = 0),
##'                 tspan = 1:10)
##'
##' ## View the C code.
##' C_code(model)
C_code <- function(model) {
    check_model_argument(model)
    model@C_code
}