Merge pull request #1149 from metrumresearchgroup/evtools

Fill out plugin for tools to implement events inside the model

R/model_include.R

@@ -174,7 +174,8 @@ plugins[["autodec"]] <- list(
 )
 
 plugins[["evtools"]] <- list(
-  name = "evtools", code = "#define _MRGSOLVE_USING_EVTOOLS_"  
+  name = "evtools", 
+  using = '#include "mrgsolve-evtools.h"'
 )
 
 # nocov end

inst/base/databox_cpp.h

@@ -46,4 +46,9 @@ double databox::tad() {
   return told < 0 ? -1.0 : time - told;
 }
 
+void databox::push(mrgsolve::evdata x) {
+  mevector.push_back(x); 
+  return;
+}
+
 #endif

inst/base/modelheader.h

@@ -98,11 +98,6 @@ typedef double capture;
 #define SQRT(a) sqrt(a)
 #endif
 
-// EVTOOLS PLUGIN
-#ifdef _MRGSOLVE_USING_EVTOOLS_
-#include "mrgsolve-evtools.h"
-#endif
-
 // These are the fundamental macros for
 // bioavailability, infusion rate, infusion duration
 // and dose lag time.  Keep these here, but

inst/base/mrgsolv.h

@@ -110,6 +110,7 @@ class databox {
   void stop_id_cf(){SYSTEMOFF=1;}///< stops solving for the current id, filling last value
   std::vector<mrgsolve::evdata> mevector;///< a collection of model events to pass back
   void mevent(double time, int evid);///< constructor for evdata objects
+  void push(mrgsolve::evdata x);
   double mtime(double time);///< creates evdata object for simple model event time
   double tad();///< calculates time after dose
 }; 

inst/base/mrgsolve-evtools.h

@@ -1,7 +1,61 @@
 
+namespace evt {
+
+typedef mrgsolve::evdata ev;
+
+mrgsolve::evdata bolus(const double amt, const int cmt) {
+  mrgsolve::evdata ev(0, 1); 
+  ev.amt = amt; 
+  ev.cmt = cmt;
+  ev.check_unique = false;
+  ev.now = true;
+  return ev;
+}
+
+void bolus(databox& self, const double amt, const int cmt) {
+  mrgsolve::evdata ev = bolus(amt, cmt);
+  self.mevector.push_back(ev);
+  return;
+}
+
+mrgsolve::evdata infuse(const double amt, const int cmt, const double rate) {
+  mrgsolve::evdata ev(0, 1); 
+  ev.amt = amt; 
+  ev.cmt = cmt;
+  ev.rate = rate;
+  ev.now = true;
+  ev.check_unique = false;
+  return ev;
+}
+
+void infuse(databox& self, const double amt, const int cmt, const double rate) {
+  mrgsolve::evdata ev = infuse(amt, cmt, rate);
+  self.mevector.push_back(ev);
+  return;
+}
+
+void retime(mrgsolve::evdata& ev, const double time) {
+  ev.time = time;
+  ev.now = false;
+  return;
+}
+
+void now(mrgsolve::evdata& ev) {
+  ev.now = true;
+  return;
+}
+
+void push(databox& self, mrgsolve::evdata ev) {
+  self.push(ev);
+  return;
+}
+
+}
+
 namespace mrgsolve {
-  namespace evt {
-    void push(databox& self, mrg::evdata ev) {
-      self.mevector.push_back(ev);  
-    }
+namespace evt {
+void push(databox& self, mrg::evdata ev) {
+  self.mevector.push_back(ev);  
+}
+}
 }

inst/maintenance/unit-cpp/test-evtools.R

@@ -0,0 +1,105 @@
+
+library(testthat)
+library(mrgsolve)
+library(dplyr)
+
+Sys.setenv(R_TESTS="")
+options("mrgsolve_mread_quiet"=TRUE)
+
+local_edition(3)
+
+code <- '
+$SET end = 12, rtol = 1e-4, atol = 1e-4, delta = 0.25
+$PLUGIN evtools
+$PARAM 
+mode = 0, f1 = 1, dose = 100, irate = 50, newtime = 2
+$CMT A B
+$PK
+F_A = f1;
+$DES
+dxdt_A = -1 * A;
+dxdt_B =  1 * A - 0.1 * B;
+$TABLE
+bool givedose = TIME==0;
+if(mode==1 && givedose) {
+  evt::bolus(self, dose, 1);
+}
+if(mode==2 && givedose) {
+  evt::infuse(self, dose, 1, irate);
+}
+if(mode==3 && givedose) {
+  evt::ev ev = evt::bolus(dose, 1); 
+  evt::retime(ev, 2);
+  evt::push(self, ev);
+}
+if(mode==4 && givedose) {
+  evt::ev ev = evt::infuse(dose, 1, irate); 
+  evt::retime(ev, newtime);
+  evt::push(self, ev);
+}
+if(mode==5 && givedose) {
+  evt::bolus(self, dose, 1); 
+  evt::bolus(self, dose, 1);
+  evt::bolus(self, dose, 1);
+}
+if(mode==6 && givedose) {
+  evt::ev ev = evt::bolus(dose, 1);
+  ev.time = 100;
+  evt::now(ev);
+  evt::push(self, ev);
+}
+'
+
+mod <- mcode("test-evtools-model-1", code)
+
+mrgsim(mod, param = list(mode = 1, f1 = 1)) 
+
+test_that("evtools - bolus now", {
+
+  out <- mrgsim(mod, param = list(mode = 1))
+  expect_equal(out$A[1], mod$dose)
+  expect_equal(out$B[1], 0)
+
+  out <- mrgsim(mod, param = list(mode = 1, f1 = 0.5))
+  expect_equal(out$A[1], mod$dose * param(out)$f1)
+  expect_equal(out$B[1], 0)
+})
+
+test_that("evtools - infusion now", {
+  out <- mrgsim(mod, param = list(mode = 2))
+  tmax <- out$time[which.max(out$A)]
+  expect_equal(tmax, mod$dose/mod$irate)
+  expect_equal(out$B[1], 0)
+
+  out <- mrgsim(mod, param = list(mode = 2, f1 = 0.5))
+  tmax <- out$time[which.max(out$A)]
+  env <- as.list(param(out)) 
+  expect_equal(tmax, param(out)$f1 * mod$dose / mod$irate)
+})
+
+test_that("evtools - bolus with retime", {
+  out <- mrgsim(mod, param = list(mode = 3), recsort = 3)
+  cmax <- filter_sims(out, time==2)
+  expect_equal(nrow(cmax), 1L)
+  expect_equal(cmax$A, mod$dose)
+})
+
+test_that("evtools - infusion with retime", {
+  out <- mrgsim(mod, param = list(mode = 4), recsort = 3)
+  tmax <- out$time[which.max(out$A)]
+  expect_equal(tmax, mod$newtime + mod$dose/mod$irate)
+  expect_equal(out$B[1], 0)
+})
+
+test_that("evtools - multiple bolus doses given now", {
+  out <- mrgsim(mod, param = list(mode = 5))
+  expect_equal(out$A[1], 300)
+  tmax <- out$time[which.max(out$A)]
+  expect_equal(tmax, 0)
+})
+
+test_that("evtools - give timed dose now", {
+  a <- mrgsim(mod, param = list(mode = 6))
+  b <- mrgsim(mod, param = list(mode = 1))
+  expect_identical(as.data.frame(a), as.data.frame(b))
+})