metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
GNU General Public License v2.0
131 stars 36 forks source link

Add $EVENT block #1230

Closed kylebaron closed 2 months ago

kylebaron commented 2 months ago

This PR adds a new block ($EVENT) which functions exactly like ($TABLE), but it gets called before $TABLE gets called. The only reason for this is so that dynamic dosing events can happen prior to output calculations in $TABLE.

Example with $EVENT

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

code <- '
$PLUGIN evtools

$CMT B

$MAIN 
double a = 1; 

$ODE
dxdt_B = -0.1*B;

$EVENT
if(TIME==0) {
  evt::ev dose = evt::infuse(100, 1, 10);
  evt::ii(dose, 24); 
  evt::addl(dose, 4);
  evt::ss(dose, 1);
  self.push(dose);
}
a = 2;

$ERROR
capture cp = B/12; 

$CAPTURE a
'

mod <- mcode("test", code, soloc = "build")
#> Creating build directory: build
#> Building test ... done.

out <- mrgsim(mod, end = 196)
plot(out)

out
#> Model:  test 
#> Dim:    197 x 5 
#> Time:   0 to 196 
#> ID:     1 
#>     ID time     B a    cp
#> 1:   1    0 17.14 2 1.429
#> 2:   1    1 25.03 2 2.086
#> 3:   1    2 32.16 2 2.680
#> 4:   1    3 38.62 2 3.218
#> 5:   1    4 44.46 2 3.705
#> 6:   1    5 49.74 2 4.145
#> 7:   1    6 54.53 2 4.544
#> 8:   1    7 58.85 2 4.905

as.list(mod)$cpp_variables
#>      type var context
#> 1  double   a    main
#> 2 capture  cp   table
p <- deparse(mod@shlib$pointers)

Example with $ERROR

code2 <- '
$PLUGIN evtools

$CMT B

$ODE
dxdt_B = -0.1*B;

$ERROR
if(TIME==0) {
  evt::ev dose = evt::infuse(100, 1, 10);
  evt::ii(dose, 24); 
  evt::addl(dose, 4);
  evt::ss(dose, 1);
  self.push(dose);
}

capture cp = B/12; 
'

mod2 <- mcode("test2", code2, soloc = "build")
#> Building test2 ... done.
out2 <- mrgsim(mod2, end = 196)
plot(out2)

out2
#> Model:  test2 
#> Dim:    197 x 4 
#> Time:   0 to 196 
#> ID:     1 
#>     ID time     B    cp
#> 1:   1    0 17.14 0.000
#> 2:   1    1 25.03 2.086
#> 3:   1    2 32.16 2.680
#> 4:   1    3 38.62 3.218
#> 5:   1    4 44.46 3.705
#> 6:   1    5 49.74 4.145
#> 7:   1    6 54.53 4.544
#> 8:   1    7 58.85 4.905

Identical outputs

head(out)
#>   ID time        B a       cp
#> 1  1    0 17.14309 2 1.428591
#> 2  1    1 25.02797 2 2.085664
#> 3  1    2 32.16250 2 2.680208
#> 4  1    3 38.61809 2 3.218174
#> 5  1    4 44.45935 2 3.704946
#> 6  1    5 49.74474 2 4.145395
head(out2)
#>   ID time        B       cp
#> 1  1    0 17.14309 0.000000
#> 2  1    1 25.02797 2.085664
#> 3  1    2 32.16250 2.680208
#> 4  1    3 38.61809 3.218174
#> 5  1    4 44.45935 3.704946
#> 6  1    5 49.74474 4.145395

identical(out$B, out2$B)
#> [1] TRUE

Created on 2024-09-10 with reprex v2.1.1

kyleam commented 2 months ago

This really isn't bad.

I agree. I think this all looks straightforward and manageable.