R/C-generator.R

In stewid/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 -- 2024 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.
##' @noRd
C_heading <- function() {
    c(sprintf("/* Generated by SimInf (v%s) */",
              utils::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",
      "")
}

C_enumeration_constants <- function(target, labels) {
    if (length(labels) == 0)
        return(character(0))

    lines <- c(
        paste0("/* Enumeration constants for indicies in the '",
               target, "' vector. */"),
        "enum {")

    for (i in seq_len(length(labels))) {
        lbl <- toupper(labels[i])
        if ((nchar(lines[length(lines)]) + nchar(lbl)) > 64)
            lines <- c(lines, "     ")
        if (i > 1)
            lines[length(lines)] <- paste0(lines[length(lines)], " ")
        lines[length(lines)] <- paste0(lines[length(lines)], lbl)
        if (i < length(labels))
            lines[length(lines)] <- paste0(lines[length(lines)], ",")
    }

    lines[length(lines)] <- paste0(lines[length(lines)], "};")
    c(lines, "")
}

C_enum <- function(compartments, ldata_names, gdata_names, v0_names,
                   use_enum) {
    if (!isTRUE(use_enum))
        return(character(0))

   c(C_enumeration_constants("u", compartments),
     C_enumeration_constants("v", v0_names),
     C_enumeration_constants("ldata", ldata_names),
     C_enumeration_constants("gdata", gdata_names))
}

##' Generate C code for the variables
##'
##' Generate C code for the variables in the model transition rate
##' functions.
##'
##' @param propensity data for the propensity.
##' @param propensity data for the variables.
##' @return character vector with C code.
##' @noRd
C_variables <- function(propensity, variables) {
    ## Determine variables to include, both directly in the
    ## propensity, and as a dependency between variables.
    j <- match(propensity$variables, names(variables))
    i <- integer(0)
    while (length(j)) {
        if (!any(i == j[1])) {
            i <- sort(c(i, j[1]))
            j <- c(j, match(variables[[j[1]]]$depends, names(variables)))
        }
        j <- j[-1]
    }

    lines <- character(0)
    for (j in i) {
        lines <- c(
            lines,
            sprintf("    const %s %s = %s;",
                    variables[[j]]$type,
                    variables[[j]]$variable,
                    variables[[j]]$code))
    }

    if (length(lines))
        lines <- c(lines, "")

    lines
}

##' 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) {
    propensities <- transitions$propensities
    variables <- transitions$variables

    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(propensities))) {
        propensity <- propensities[[i]]

        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,
            "{")

        ## Insert variables.
        lines <- c(
            lines,
            C_variables(propensity, variables))

        ## Insert propensity.
        lines <- c(
            lines,
            sprintf("    return %s;", propensity$code),
            "}",
            "")
    }

    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) {
    n_tr_fn <- length(transitions$propensities)

    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(n_tr_fn), 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, compartments,
                          ldata_names, gdata_names, v0_names,
                          use_enum) {
    c(C_heading(),
      C_include(),
      C_define(),
      C_enum(compartments, ldata_names, gdata_names, v0_names,
             use_enum),
      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
}