Providing Custom C/C++ Functions to rxode2 from a Package

Overview

This vignette explains how to write an R package that provides custom C or C++ functions for use inside rxode2 (and nlmixr2) models. Compared to the inline approach, packaging the C code lets you:

provide proper documentation, tests, and versioning for the functions that are not yet provided in nlmixr2/rxode2.
supply exact analytic derivatives so that nlmixr2 optimization works without falling back to finite differences.

The worked example uses a package named mm that implements a simple Michaelis-Menten saturable-elimination function

$\text{mm}(C, V_{\max}, K_m) = \frac{V_{\max} \cdot C}{K_m + C}$

along with its first- and second-order partial derivatives.

How rxode2 calls compiled functions

When rxode2 compiles a model it generates a C source file. Inside that file, every external function the model uses is resolved at load time with a call like

mm = (rxode2_fn3) R_GetCCallable("mm", "mm");

R_GetCCallable is R’s standard mechanism for C-level inter-package calls. For it to succeed, your package must register the symbol with R_RegisterCCallable during its own initialization.

The mapping from a model expression (mm(Cp, Vmax, Km)) to the generated R_GetCCallable call is held in a data frame that rxode2 maintains at run-time. You add rows to that data frame in .onAttach using rxode2::rxode2parseAssignTranslation().

Package skeleton

mm/
├── DESCRIPTION
├── NAMESPACE
├── R/
│   ├── mm.R           # R interface + roxygen docs
│   └── zzz.R          # .onAttach / .onUnload hooks
└── src/
    ├── myfuns.c        # scalar C functions and SEXP vector wrappers
    └── init.c          # R_init_mm: registers callables and routines

DESCRIPTION

Package: mm
Version: 0.1.0
Imports: rxode2
LinkingTo: rxode2
RoxygenNote: 7.3.2

LinkingTo: rxode2 gives your C code access to the rxode2 headers, including the rxode2_fn* typedef family used in the translation table. RoxygenNote is required when using roxygen2 to generate documentation.

NAMESPACE

The NAMESPACE is generated by roxygen2 from tags in R/mm.R and R/zzz.R; you do not edit it by hand. The relevant tags are shown with each file below.

Writing the C functions

`src/myfuns.c`

The file contains two layers:

Scalar functions (double in, double out) — these are what rxode2 calls inside the ODE solver.
SEXP vector wrappers (SEXP in, SEXP out) — these are what the thin R interface calls; they loop over equal-length double vectors and return the result as an R numeric vector.

/* Standard R headers ---------------------------------------- */
#define STRICT_R_HEADERS
#include <R.h>
#include <Rinternals.h>
#include <R_ext/RS.h>

/* ============================================================
   Scalar functions (used by rxode2 ODE solver)
   ============================================================ */

/* mm(C, Vmax, Km) = Vmax*C / (Km+C) */
double mm(double C, double Vmax, double Km) {
  return Vmax * C / (Km + C);
}

/* d(mm)/dC   = Vmax*Km / (Km+C)^2 */
double mm_dC(double C, double Vmax, double Km) {
  double den = Km + C;
  return Vmax * Km / (den * den);
}

/* d(mm)/dVmax = C / (Km+C) */
double mm_dVmax(double C, double Vmax, double Km) {
  (void)Vmax;
  return C / (Km + C);
}

/* d(mm)/dKm  = -Vmax*C / (Km+C)^2 */
double mm_dKm(double C, double Vmax, double Km) {
  double den = Km + C;
  return -Vmax * C / (den * den);
}

/* ============================================================
   SEXP vector wrappers (called from R via .Call)
   Assumes C, Vmax, Km are REALSXP vectors of the same length
   (enforce equal length via data.frame() in the R interface).
   ============================================================ */

SEXP _mm_mm(SEXP C, SEXP Vmax, SEXP Km) {
  int n = LENGTH(C);
  SEXP out = PROTECT(allocVector(REALSXP, n));
  double *c = REAL(C), *vmax = REAL(Vmax), *km = REAL(Km);
  double *res = REAL(out);
  for (int i = 0; i < n; i++) res[i] = mm(c[i], vmax[i], km[i]);
  UNPROTECT(1);
  return out;
}

SEXP _mm_mm_dC(SEXP C, SEXP Vmax, SEXP Km) {
  int n = LENGTH(C);
  SEXP out = PROTECT(allocVector(REALSXP, n));
  double *c = REAL(C), *vmax = REAL(Vmax), *km = REAL(Km);
  double *res = REAL(out);
  for (int i = 0; i < n; i++) res[i] = mm_dC(c[i], vmax[i], km[i]);
  UNPROTECT(1);
  return out;
}

SEXP _mm_mm_dVmax(SEXP C, SEXP Vmax, SEXP Km) {
  int n = LENGTH(C);
  SEXP out = PROTECT(allocVector(REALSXP, n));
  double *c = REAL(C), *vmax = REAL(Vmax), *km = REAL(Km);
  double *res = REAL(out);
  for (int i = 0; i < n; i++) res[i] = mm_dVmax(c[i], vmax[i], km[i]);
  UNPROTECT(1);
  return out;
}

SEXP _mm_mm_dKm(SEXP C, SEXP Vmax, SEXP Km) {
  int n = LENGTH(C);
  SEXP out = PROTECT(allocVector(REALSXP, n));
  double *c = REAL(C), *vmax = REAL(Vmax), *km = REAL(Km);
  double *res = REAL(out);
  for (int i = 0; i < n; i++) res[i] = mm_dKm(c[i], vmax[i], km[i]);
  UNPROTECT(1);
  return out;
}

Rules every rxode2-compatible scalar C function must follow:

Return type is double.
All arguments are double scalars, not vectors.
No R API calls inside the body — pure numeric C.

`src/init.c` — registering the callables

#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>
#include <stdlib.h>   /* NULL */

/* Forward-declare scalar functions from myfuns.c */
extern double mm      (double, double, double);
extern double mm_dC   (double, double, double);
extern double mm_dVmax(double, double, double);
extern double mm_dKm  (double, double, double);

/* Forward-declare SEXP wrappers from myfuns.c */
extern SEXP _mm_mm      (SEXP, SEXP, SEXP);
extern SEXP _mm_mm_dC   (SEXP, SEXP, SEXP);
extern SEXP _mm_mm_dVmax(SEXP, SEXP, SEXP);
extern SEXP _mm_mm_dKm  (SEXP, SEXP, SEXP);

void R_init_mm(DllInfo *dll) {

  /* --- .Call registration -----------------------------------
     Makes SEXP wrappers available to .Call() from R.          */
  static const R_CallMethodDef callMethods[] = {
    {"_mm_mm",       (DL_FUNC) &_mm_mm,       3},
    {"_mm_mm_dC",    (DL_FUNC) &_mm_mm_dC,    3},
    {"_mm_mm_dVmax", (DL_FUNC) &_mm_mm_dVmax, 3},
    {"_mm_mm_dKm",   (DL_FUNC) &_mm_mm_dKm,  3},
    {NULL, NULL, 0}
  };
  R_registerRoutines(dll, NULL, callMethods, NULL, NULL);
  R_useDynamicSymbols(dll, FALSE);

  /* --- R_RegisterCCallable ----------------------------------
     REQUIRED for rxode2 integration.  Makes each scalar symbol
     retrievable from other packages' C code via
     R_GetCCallable("mm", "<name>").
     rxode2's generated model C files use exactly this.        */
  R_RegisterCCallable("mm", "mm",       (DL_FUNC) &mm);
  R_RegisterCCallable("mm", "mm_dC",    (DL_FUNC) &mm_dC);
  R_RegisterCCallable("mm", "mm_dVmax", (DL_FUNC) &mm_dVmax);
  R_RegisterCCallable("mm", "mm_dKm",   (DL_FUNC) &mm_dKm);
}

R_RegisterCCallable and R_registerRoutines serve different purposes:

Function	Who uses it	Purpose
`R_registerRoutines`	R itself	`.Call()` from R code
`R_RegisterCCallable`	Other packages’ compiled C	`R_GetCCallable(pkg, sym)`

rxode2 uses R_GetCCallable, so R_RegisterCCallable is the critical call for rxode2 integration.

R interface

`R/mm.R`

The R interface is optional but useful: it lets analysts call mm() directly on vectors (e.g. in dplyr pipelines or plots) without needing to understand the C internals. The data.frame() call handles recycling and length checking automatically.

#' Michaelis-Menten saturable elimination
#'
#' Computes \eqn{V_{\max} \cdot C / (K_m + C)} element-wise over
#' equal-length (or recyclable) numeric vectors.
#'
#' @param C Numeric vector of concentrations.
#' @param Vmax Numeric vector of maximum elimination rates.
#' @param Km Numeric vector of Michaelis constants.
#' @return Numeric vector of elimination rates, the same length as the
#'   recycled inputs.
#' @export
#' @useDynLib mm, .registration = TRUE
#' @importFrom rxode2 rxD rxode2parseAssignTranslation rxode2parseGetTranslation
mm <- function(C, Vmax, Km) {
  df <- data.frame(C = C, Vmax = Vmax, Km = Km)
  .Call(`_mm_mm`, as.double(df$C), as.double(df$Vmax), as.double(df$Km))
}

#' @describeIn mm Partial derivative with respect to C
#' @export
mm_dC <- function(C, Vmax, Km) {
  df <- data.frame(C = C, Vmax = Vmax, Km = Km)
  .Call(`_mm_mm_dC`, as.double(df$C), as.double(df$Vmax), as.double(df$Km))
}

#' @describeIn mm Partial derivative with respect to Vmax
#' @export
mm_dVmax <- function(C, Vmax, Km) {
  df <- data.frame(C = C, Vmax = Vmax, Km = Km)
  .Call(`_mm_mm_dVmax`, as.double(df$C), as.double(df$Vmax), as.double(df$Km))
}

#' @describeIn mm Partial derivative with respect to Km
#' @export
mm_dKm <- function(C, Vmax, Km) {
  df <- data.frame(C = C, Vmax = Vmax, Km = Km)
  .Call(`_mm_mm_dKm`, as.double(df$C), as.double(df$Vmax), as.double(df$Km))
}

Registering with rxode2

`R/zzz.R`

.onAttach <- function(libname, pkgname) {

  ## 1. Add rows to rxode2's function translation table.
  ##
  ##    Column meanings:
  ##      rxFun       - name used in rxode2 model code
  ##      fun         - C variable name in the generated model source
  ##      type        - rxode2 function-pointer typedef (call convention)
  ##      package     - first argument to R_GetCCallable()
  ##      packageFun  - second argument to R_GetCCallable()
  ##      argMin      - minimum number of accepted arguments
  ##      argMax      - maximum number of accepted arguments
  ##      threadSafe  - 1 if safe to call from OpenMP threads, 0 otherwise
  ##
  ##    Note: argMin and argMax can only differ when the function parsing
  ##    is directly handled by rxode2 internals.  For package-provided
  ##    functions they must be equal.
  ##
  ##    Common function-pointer typedefs (rxode2_model_shared.h):
  ##      rxode2_fn    double (*)(double)
  ##      rxode2_fn2   double (*)(double, double)
  ##      rxode2_fn3   double (*)(double, double, double)
  ##
  ##    All four functions here take three double args -> "rxode2_fn3".
  ##    Registering the derivative functions in the table allows rxode2
  ##    to call them directly as C when assembling the Jacobian.
  ##
  ##    Only add the rows if this package's functions are not already
  ##    present (e.g. after a package reload in an interactive session).

  .cur <- rxode2::rxode2parseGetTranslation()

  .newRows <- data.frame(
    rxFun      = c("mm",         "mm_dC",      "mm_dVmax",   "mm_dKm"),
    fun        = c("mm",         "mm_dC",      "mm_dVmax",   "mm_dKm"),
    type       = rep("rxode2_fn3", 4),
    package    = rep("mm", 4),
    packageFun = c("mm",         "mm_dC",      "mm_dVmax",   "mm_dKm"),
    argMin     = 3L,
    argMax     = 3L,
    threadSafe = 1L,
    stringsAsFactors = FALSE
  )

  ## Check for name conflicts before modifying the table.
  .mmNames <- .newRows$rxFun
  .conflicts <- .cur[.cur$rxFun %in% .mmNames, , drop = FALSE]

  if (nrow(.conflicts) > 0L) {
    .fromOther <- .conflicts[.conflicts$package != "mm", , drop = FALSE]
    .fromSelf  <- .conflicts[.conflicts$package == "mm", , drop = FALSE]

    if (nrow(.fromOther) > 0L) {
      stop(
        "Cannot load package 'mm': the following rxode2 function name(s) are ",
        "already registered by another package:\n",
        paste0("  ", .fromOther$rxFun, " (registered by '",
               .fromOther$package, "')", collapse = "\n"),
        call. = FALSE
      )
    }

    if (nrow(.fromSelf) > 0L) {
      packageStartupMessage(
        "Package 'mm' rxode2 function registration was already present ",
        "(package reloaded); re-registering."
      )
    }
  }

  ## Remove any existing mm rows (covers the reload case) then add fresh ones.
  .cur <- .cur[.cur$package != "mm", , drop = FALSE]

  rxode2::rxode2parseAssignTranslation(rbind(.cur, .newRows))

  ## 2. Register derivatives for nlmixr2 gradient computation.
  ##
  ##    rxD(name, list-of-derivative-functions)
  ##
  ##    Each function receives the same argument names as the original
  ##    and must return a CHARACTER string: the C expression for the
  ##    partial derivative w.r.t. that argument.
  ##
  ##    For mm() we reference the registered C derivative functions
  ##    directly rather than writing out the algebra inline.  This
  ##    keeps the expressions short and delegates correctness to the
  ##    C implementation.

  rxode2::rxD("mm", list(
    function(C, Vmax, Km) paste0("mm_dC(",    C, ",", Vmax, ",", Km, ")"),
    function(C, Vmax, Km) paste0("mm_dVmax(", C, ",", Vmax, ",", Km, ")"),
    function(C, Vmax, Km) paste0("mm_dKm(",   C, ",", Vmax, ",", Km, ")")
  ))

  ## 3. Register second-order derivatives as inline C expressions.
  ##    These are the derivatives of mm_dC, mm_dVmax, and mm_dKm,
  ##    supplied as expression strings rather than additional C functions.
  ##
  ##    mm_dC  = Vmax*Km/(Km+C)^2
  ##      d/dC    = -2*Vmax*Km / (Km+C)^3
  ##      d/dVmax =    Km      / (Km+C)^2
  ##      d/dKm   =  Vmax*(C-Km) / (Km+C)^3

  rxode2::rxD("mm_dC", list(
    function(C, Vmax, Km)
      paste0("(-2.0*(", Vmax, ")*(", Km, "))/pow((", Km, ")+(", C, "),3.0)"),
    function(C, Vmax, Km)
      paste0("(", Km, ")/pow((", Km, ")+(", C, "),2.0)"),
    function(C, Vmax, Km)
      paste0("(", Vmax, ")*((", C, ")-(", Km, "))/pow((", Km, ")+(", C, "),3.0)")
  ))

  ##    mm_dVmax = C/(Km+C)
  ##      d/dC    =  Km / (Km+C)^2
  ##      d/dVmax =  0
  ##      d/dKm   = -C  / (Km+C)^2

  rxode2::rxD("mm_dVmax", list(
    function(C, Vmax, Km)
      paste0("(", Km, ")/pow((", Km, ")+(", C, "),2.0)"),
    function(C, Vmax, Km) "0.0",
    function(C, Vmax, Km)
      paste0("-(", C, ")/pow((", Km, ")+(", C, "),2.0)")
  ))

  ##    mm_dKm  = -Vmax*C/(Km+C)^2
  ##      d/dC    =  Vmax*(C-Km) / (Km+C)^3
  ##      d/dVmax = -C           / (Km+C)^2
  ##      d/dKm   =  2*Vmax*C   / (Km+C)^3

  rxode2::rxD("mm_dKm", list(
    function(C, Vmax, Km)
      paste0("(", Vmax, ")*((", C, ")-(", Km, "))/pow((", Km, ")+(", C, "),3.0)"),
    function(C, Vmax, Km)
      paste0("-(", C, ")/pow((", Km, ")+(", C, "),2.0)"),
    function(C, Vmax, Km)
      paste0("2.0*(", Vmax, ")*(", C, ")/pow((", Km, ")+(", C, "),3.0)")
  ))

  invisible()
}

.onUnload <- function(libpath) {
  ## Remove only this package's rows from the translation table.
  ## Filtering by package == "mm" is safe when multiple packages have
  ## each added their own rows: other packages' entries are untouched.
  .cur <- try(rxode2::rxode2parseGetTranslation(), silent = TRUE)
  if (!inherits(.cur, "try-error")) {
    try(rxode2::rxode2parseAssignTranslation(
      .cur[.cur$package != "mm", , drop = FALSE]
    ), silent = TRUE)
  }
  ## Remove derivative entries for all registered function names.
  try(rxode2::rxRmFun("mm"),       silent = TRUE)
  try(rxode2::rxRmFun("mm_dC"),    silent = TRUE)
  try(rxode2::rxRmFun("mm_dVmax"), silent = TRUE)
  try(rxode2::rxRmFun("mm_dKm"),   silent = TRUE)
}

How it works end-to-end

When a user writes rxode2({ dA1/dt <- -mm(Cp, Vmax, Km) }):

rxode2 looks up "mm" in the translation data frame and finds type = "rxode2_fn3", package = "mm", packageFun = "mm".

It generates C code containing:

mm       = (rxode2_fn3) R_GetCCallable("mm", "mm");
mm_dC    = (rxode2_fn3) R_GetCCallable("mm", "mm_dC");
mm_dVmax = (rxode2_fn3) R_GetCCallable("mm", "mm_dVmax");
mm_dKm   = (rxode2_fn3) R_GetCCallable("mm", "mm_dKm");

When the compiled model DLL is loaded these calls resolve to the functions exported from the mm package via R_RegisterCCallable.
Calls to mm() inside the ODE system execute directly as C with no R overhead and full OpenMP thread safety.
When nlmixr2 needs the gradient it calls the derivative expressions registered by rxD("mm", ...), which expand to mm_dC(...), mm_dVmax(...), and mm_dKm(...) — also resolved directly as C.

Named and variable-length arguments via `rxUdfUi`

The C interface requires arguments in a fixed positional order (C, Vmax, Km). rxode2 models are translated to C, so named arguments and alternative argument orders are not supported at the C level. However, you can accept them at the model-parsing level by providing an rxUdfUi.mm S3 method that rewrites the call before it reaches the C translator.

This mechanism is described in full in the Integrating User Defined Functions into rxode2 vignette (section Functions to insert rxode2 code into the current model). The short version: rxode2 calls rxUdfUi.mm(fun) during model parsing whenever it encounters a call to mm() in a ui model, where fun is the quoted language object of that call. The method must return a list with at minimum a $replace element containing the rewritten call as a character string in the canonical positional form that the C translator expects.

`R/mm.R` — adding `rxUdfUi.mm`


#' @importFrom rxode2 rxUdfUi
#' @export
rxode2::rxUdfUi

#' rxode2 ui model interface for mm()
#'
#' Allows \code{mm()} to be called with named arguments in any order
#' inside \code{rxode2}/\code{nlmixr2} ui models, e.g.
#' \code{mm(Km=Km, C=Cp, Vmax=Vmax)}.  The method rewrites the call
#' to the canonical positional form \code{mm(C, Vmax, Km)} before the
#' model is compiled.
#'
#' @param fun Quoted language object representing the \code{mm()} call
#'   as it appears in the model.
#' @return A list with element \code{$replace}: the rewritten call as
#'   a character string.
#' @export
rxUdfUi.mm <- function(fun) {
  # Use a dummy function with the canonical argument order so that
  # match.call() can map any combination of named / positional
  # arguments back to the C-level positional order.
  .dummy <- function(C, Vmax, Km) {}
  .mc <- match.call(.dummy, fun)

  .C    <- deparse1(.mc$C)
  .Vmax <- deparse1(.mc$Vmax)
  .Km   <- deparse1(.mc$Km)

  list(replace = paste0("mm(", .C, ",", .Vmax, ",", .Km, ")"))
}

Exporting rxUdfUi.mm with @export is sufficient for roxygen2 to register the S3 method when the package is installed — no .s3register call is needed in .onAttach.

With this method in place, all of the following are equivalent inside a ui model:

## Positional (always worked):
m1 <- function() {
  model({
    d/dt(A1) <- -mm(Cp, Vmax, Km)
  })
}


## Named in canonical order:
m2 <- function() {
  model({
    d/dt(A1) <- -mm(C = Cp, Vmax = Vmax, Km = Km)
  })
}

## Named in a different order — rxUdfUi.mm reorders before compilation:
m3 <- function() {
  model({
    d/dt(A1) <- -mm(Km = Km, Vmax = Vmax, C = Cp)
  })
}

Note that plain rxode2({...}) (low-level, without ui) does not call rxUdfUi methods and still requires positional arguments.

Using the package

library(rxode2)
library(mm)   # triggers .onAttach -> rxode2parseAssignTranslation + rxD

# rxode2 produces a syntax error for the wrong number of arguments:
try(rxode2({ d/dt(A1) <- -mm(Cp, Vmax) }))


## --- rxode2 simulation ---
mod <- rxode2({
  Cp       <- A1 / Vc
  d/dt(A1) <- -mm(Cp, Vmax, Km)
})

ev  <- et(amt = 100) |> et(0, 24, by = 0.5)
sim <- rxSolve(mod, ev, params = c(Vc = 10, Vmax = 5, Km = 1))
plot(sim, Cp)

## --- R vector interface ---
mm(C = c(1, 2, 5), Vmax = 10, Km = 2)   # vectorised over C

Verifying registration

These could be adapted to tests so that you can make sure that rxode2 loads the external function as requested:

## Check translation table
trans <- rxode2::rxode2parseGetTranslation()
trans[trans$rxFun %in% c("mm", "mm_dC", "mm_dVmax", "mm_dKm"), ]

## Check first-order derivatives (should return C function-call strings)
rxode2::rxFromSE("Derivative(mm(C, Vmax, Km), C)")
rxode2::rxFromSE("Derivative(mm(C, Vmax, Km), Vmax)")
rxode2::rxFromSE("Derivative(mm(C, Vmax, Km), Km)")

## Check second-order derivatives (should return inline C expressions)
rxode2::rxFromSE("Derivative(mm_dC(C, Vmax, Km), C)")
rxode2::rxFromSE("Derivative(mm_dKm(C, Vmax, Km), Km)")

Summary checklist

Step	File	What to do
1	`src/myfuns.c`	Write scalar `double fun(double, ...)` functions
2	`src/myfuns.c`	Write `SEXP _pkg_fun(SEXP, ...)` vector wrappers calling the scalars
3	`src/init.c`	Register SEXP wrappers with `R_registerRoutines` (`.Call` table)
4	`src/init.c`	Register scalar functions with `R_RegisterCCallable`
5	`DESCRIPTION`	Add `Imports: rxode2`, `LinkingTo: rxode2`, `RoxygenNote:`
6	`R/mm.R`	Write thin R wrappers using `.Call(\`_pkg_fun`, …)`with roxygen docs \| \| 6b \|`R/mm.R`\| Optionally add`rxUdfUi.mm`to support named / reordered arguments in ui models \| \| 7 \|`R/zzz.R`\| Call`rxode2parseAssignTranslation(rbind(old, newRows))`in`.onAttach`\| \| 8 \|`R/zzz.R`\| Call`rxD(“fun”, list(…))`referencing C functions for 1st-order \| \| 9 \|`R/zzz.R`\| Call`rxD(“fun_dArg”, list(…))`with expression strings for 2nd-order \| \| 10 \|`R/zzz.R`\| Filter out rows where`package == “mm”`+ call`rxRmFun`for all names in`.onUnload`