728 wnl style linear compartment models using combined rxode2

NAMESPACE

@@ -286,6 +286,8 @@ export(.rxWithSink)
 export(.rxWithSinkBoth)
 export(.rxWithWd)
 export(.s3register)
+export(.solComp2)
+export(.solComp3)
 export(.symengineFs)
 export(.toClassicEvid)
 export(.udfCallFunArg)

R/linCmt.R

@@ -0,0 +1,56 @@
+#' Calculate the lambdas and coefficients of the two compartment model
+#'
+#' @param k10 elimination rate
+#' @param k12 rate from central to peripheral compartment
+#' @param k21 rate from peripheral to central compartment
+#' @param cpp boolean; call the c++ interface used for stan gradients
+#' @return List with `L` vector and matrices `C1` and `C2`
+#' @export
+#' @keywords internal
+#' @author Matthew L. Fidler based on `wnl` package/paper, implemented
+#'   in C/C++
+#' @examples
+#' .solComp2(k10=0.1, k12=3, k21=1)
+.solComp2 <- function(k10, k12, k21, cpp=FALSE) {
+  checkmate::assertNumeric(k10, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k12, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k21, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertLogical(cpp, len=1, any.missing = FALSE)
+  if (cpp) {
+    .ret <- .Call(`_rxode2_solComp2cpp`, k10, k12, k21)
+  } else {
+    .ret <- .Call(`_rxode2_solComp2`, k10, k12, k21)
+  }
+  if (is.null(.ret)) {
+    stop("roots must be distinct real values", call.=FALSE)
+  }
+  .ret
+}
+#' Calculate the lambdas and coefficients of the three compartment model
+#'
+#' @inheritParams .solComp2
+#' @param k13 rate from central to peripheral compartment #2
+#' @param k31 rate from peripheral compartment #2 to central
+#' @return List with `L` vector and matrices `C1`, `C2` and `C3`
+#' @export
+#' @author Matthew L. Fidler
+#' @keywords internal
+#' @examples
+#' .solComp3(k10=0.1, k12=3, k21=1, k13=2, k31=0.5)
+.solComp3 <- function(k10, k12, k21, k13, k31, cpp=FALSE) {
+  checkmate::assertNumeric(k10, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k12, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k21, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k13, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(k31, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertLogical(cpp, len=1, any.missing = FALSE)
+  if (cpp) {
+    .ret <- .Call(`_rxode2_solComp3cpp`, k10, k12, k21, k13, k31)
+  } else {
+    .ret <- .Call(`_rxode2_solComp3`, k10, k12, k21, k13, k31)
+  }
+  if (is.null(.ret)) {
+    stop("roots must be distinct real values", call.=FALSE)
+  }
+  .ret
+}

inst/include/rxode2.h

@@ -30,11 +30,14 @@
 #define getSolve(idx) ind->solve + (op->neq + op->nlin)*(idx)
 
 #ifdef _isrxode2_
-
+#define min2( a , b )  ( (a) < (b) ? (a) : (b) )
 #define max2( a , b )  ( (a) > (b) ? (a) : (b) )
 #define isSameTime(xout, xp) (fabs((xout)-(xp))  <= DBL_EPSILON*max2(fabs(xout),fabs(xp)))
+
 // use ~dop853 definition of same time
 #define isSameTimeDop(xout, xp) (0.1 * fabs((xout)-(xp)) <= fabs(xout) * 2.3E-16)
+#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp))
+
 
 #else
 

src/approx.c

@@ -13,8 +13,6 @@
 #define _as_zero(a) (fabs(a) < sqrt(DBL_EPSILON) ? 0.0 : a)
 #define _as_dbleps(a) (fabs(a) < sqrt(DBL_EPSILON) ? ((a) < 0 ? -sqrt(DBL_EPSILON)  : sqrt(DBL_EPSILON)) : a)
 
-#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp))
-
 #ifdef ENABLE_NLS
 #include <libintl.h>
 #define _(String) dgettext ("rxode2", String)

src/comp.c

@@ -0,0 +1,71 @@
+#define USE_FC_LEN_T
+#define STRICT_R_HEADERS
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>   /* dj: import intptr_t */
+//#include "ode.h"
+#define R_NO_REMAP
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
+#include <Rmath.h>
+#include <R_ext/Lapack.h>
+#include <R_ext/BLAS.h>
+#ifndef FCONE
+# define FCONE
+#endif
+
+#include "../inst/include/rxode2.h"
+#include "comp.h"
+#include "parTrans.h"
+
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("rxode2", String)
+/* replace pkg as appropriate */
+#else
+#define _(String) (String)
+#endif
+
+int comp1c(int ncmt,
+           int trans,
+           double *yp, // prior y
+           double xp, // initial time point
+           double xout, // final time point
+           double *p,
+           double *rate) {
+  // no update for same time
+  if (isSameTime(xout, xp)) return 1;
+  lin_context_c_t lin;
+  lin.ncmt = ncmt;
+  lin.ka=lin.k10=
+    lin.k12=lin.k21=
+    lin.k13=lin.k31=lin.v=
+    lin.dt=xout-xp;
+  lin.rate = rate;
+  if (!parTrans(&trans, &p[0], &p[1], &p[2], &p[3], &p[4], &p[5],
+                &ncmt, &(lin.k10), &(lin.v), &(lin.k12),
+                &(lin.k21), &(lin.k13), &(lin.k31))){
+    return 0;
+  }
+  return comp1(yp, &lin);
+}
+
+// This is a R wrapper for comp1c
+SEXP _rxode2_comp1c(SEXP prior, SEXP xpS, SEXP xoutS,
+                    SEXP rateS, SEXP transS, SEXP parS) {
+  SEXP retS = PROTECT(Rf_allocVector(INTSXP, 1));
+  int *ret = INTEGER(retS);
+  ret[0] = 0;
+  int ncmt = INTEGER(VECTOR_ELT(transS, 1))[0];
+  int trans = INTEGER(VECTOR_ELT(transS, 2))[0];
+  comp1c(ncmt, trans, REAL(prior), // prior y
+         REAL(xpS)[0], // initial time point
+         REAL(xoutS)[0], // final time point
+         REAL(rateS),
+         REAL(rateS));
+  UNPROTECT(1);
+  return retS;
+}