Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add $EVENT block #1230

Merged
merged 28 commits into from
Sep 25, 2024
Merged

Add $EVENT block #1230

merged 28 commits into from
Sep 25, 2024

Conversation

kylebaron
Copy link
Collaborator

@kylebaron kylebaron commented Sep 8, 2024

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

if(call_event && !this_rec->is_lagged()) {
prob.table_call();
}

if(this_rec->output()) {
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Second bit of "workaround" code to maintain event creation in both $EVENT and $TABLE.

@@ -461,6 +462,7 @@ mread <- function(model, project = getOption("mrgsolve.project", getwd()),
x@shlib[["nm_import"]] <- mread.env[["nm_import"]]
x@shlib[["source"]] <- file.path(build[["soloc"]],build[["compfile"]])
x@shlib[["md5"]] <- build[["md5"]]
x@shlib[["call_event"]] <- "EVENT" %in% names(spec)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We will pass this into the C+ simulation code to tell if we are using $EVENT to create modeled doses.

prob.event_call();
} else {
prob.table_call();
}
}
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the first work around to maintain event code in $EVENT and $TABLE. Basically, if the user has $EVENT block, then call that one here; otherwise call $TABLE. Down below, we'll call $TABLE if $EVENT was called and do nothing if it wasn't.

// MODELED EVENTS:
__BEGIN_event__
__END_event__

// TABLE CODE BLOCK:
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was programmatically generaged; every model will have this in the source going fowward.

//! signature for <code>$EVENT</code>
#define MRGSOLVE_EVENT_SIGNATURE const dvec& _A_, const dvec& _A_0_, dvec& _THETA_, const dvec& _F_, const dvec& _R_, databox& self, const dvec& _pred_, dvec& _capture_, mrgsolve::resim& simeps
#define MRGSOLVE_EVENT_SIGNATURE_N 9

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is important; I'm having $EVENT function identical to $TABLE. Anything you could do in one, you could do in the other.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could re-use that signature macro to emphasize this (e.g. use MRGSOLVE_TABLE_SIGNATURE for the event; what do you think? maybe it would be confusing or look like a mistake or it could ensure they act the same.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it could ensure they act the same.

I like the idea of sharing for this reason, and I think the confusion could be made less likely by the combination of 1) putting in a code comment and 2) moving the macros next to each other.

Another idea: what about renaming MRGSOLVE_TABLE_SIGNATURE/MRGSOLVE_TABLE_SIGNATURE_N to indicate both TABLE and EVENT and then using it for both? I suppose MRGSOLVE_TABLE_OR_EVENT_SIGNATURE is getting a bit long, but I think it'd resolve both points above.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I moved the macros next to each other and noted they are the same. I ended up copying the TABLE macro into the EVENT macro. So they will have different names in the code, but guaranteed to be in sync.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm going to send this off to winbuilder just to make sure the macros are behaving as I think they will on gcc. If that comes back ok, will merge.

@kylebaron kylebaron marked this pull request as draft September 10, 2024 14:06
@kyleam
Copy link
Contributor

kyleam commented Sep 11, 2024

This really isn't bad.

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

@@ -794,7 +794,8 @@ parin <- function(x) {
ss_n = 500,
ss_fixed = FALSE,
interrupt = -1,
etasrc = "omega"
etasrc = "omega",
call_event = x@shlib$call_event
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We've saved a logical value indicating of $EVENT is in the model; if so, tell whether or not to call.

model = "",
omats,
smats,
set = list(),
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just light touchup; this code is really old and probably don't need to pass all the functions separately.

@@ -368,7 +368,7 @@ param_re_find <- "\\bparam\\s+\\w+\\s*="
# please-deprecate
move_global <- function(x,env) {

what <- intersect(c("PREAMBLE","MAIN", "ODE", "TABLE", "PRED"),names(x))
what <- intersect(c("PREAMBLE","MAIN", "ODE", "TABLE", "EVENT", "PRED"),names(x))
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I updated this, but it's not in use any more (see comment). I'd rather get rid of it at this point, but it may or may not be in scope of this PR.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Taking a quick look at move_global, I think a dedicated/separate PR would be better.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree.

to_ns <- bind_rows(
pream$vars,
pred$vars,
main$vars,
ode$vars,
table$vars
table$vars,
event$vars
)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All of this is happening in the new version (move_global2()).

init_fun = main_func(x),
table_fun = table_func(x),
event_fun = event_func(x),
config_fun = config_func(x),
model = model(x),
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added argument names for safety ... was getting called by position before.

//! signature for <code>$EVENT</code>
#define MRGSOLVE_EVENT_SIGNATURE const dvec& _A_, const dvec& _A_0_, dvec& _THETA_, const dvec& _F_, const dvec& _R_, databox& self, const dvec& _pred_, dvec& _capture_, mrgsolve::resim& simeps
#define MRGSOLVE_EVENT_SIGNATURE_N 9

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could re-use that signature macro to emphasize this (e.g. use MRGSOLVE_TABLE_SIGNATURE for the event; what do you think? maybe it would be confusing or look like a mistake or it could ensure they act the same.

@kylebaron kylebaron marked this pull request as ready for review September 24, 2024 13:03
Copy link
Contributor

@kyleam kyleam left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM

//! signature for <code>$EVENT</code>
#define MRGSOLVE_EVENT_SIGNATURE const dvec& _A_, const dvec& _A_0_, dvec& _THETA_, const dvec& _F_, const dvec& _R_, databox& self, const dvec& _pred_, dvec& _capture_, mrgsolve::resim& simeps
#define MRGSOLVE_EVENT_SIGNATURE_N 9

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it could ensure they act the same.

I like the idea of sharing for this reason, and I think the confusion could be made less likely by the combination of 1) putting in a code comment and 2) moving the macros next to each other.

Another idea: what about renaming MRGSOLVE_TABLE_SIGNATURE/MRGSOLVE_TABLE_SIGNATURE_N to indicate both TABLE and EVENT and then using it for both? I suppose MRGSOLVE_TABLE_OR_EVENT_SIGNATURE is getting a bit long, but I think it'd resolve both points above.

Comment on lines +77 to +78
# Known that cp isn't calculated when $TABLE is used
expect_equal(outer$cp[1], 0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice.

@@ -368,7 +368,7 @@ param_re_find <- "\\bparam\\s+\\w+\\s*="
# please-deprecate
move_global <- function(x,env) {

what <- intersect(c("PREAMBLE","MAIN", "ODE", "TABLE", "PRED"),names(x))
what <- intersect(c("PREAMBLE","MAIN", "ODE", "TABLE", "EVENT", "PRED"),names(x))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Taking a quick look at move_global, I think a dedicated/separate PR would be better.

//! signature for <code>$EVENT</code> same as what we use for <code>$TABLE</code>
#define MRGSOLVE_EVENT_SIGNATURE MRGSOLVE_TABLE_SIGNATURE
#define MRGSOLVE_EVENT_SIGNATURE_N MRGSOLVE_TABLE_SIGNATURE_N

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Re ordered so that EVENT is by TABLE and re-using the definition so that it's in sync.

* using log directory 'd:/RCompile/CRANguest/R-release/mrgsolve.Rcheck'
* using R version 4.4.1 (2024-06-14 ucrt)
* using platform: x86_64-w64-mingw32
* R was compiled by
    gcc.exe (GCC) 13.2.0
    GNU Fortran (GCC) 13.2.0
* running under: Windows Server 2022 x64 (build 20348)
* using session charset: UTF-8
* checking for file 'mrgsolve/DESCRIPTION' ... OK
* checking extension type ... Package
* this is package 'mrgsolve' version '1.5.1.9002'
* package encoding: UTF-8
* checking CRAN incoming feasibility ... [11s] NOTE
Maintainer: 'Kyle T Baron <kyleb@metrumrg.com>'

Version contains large components (1.5.1.9002)
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking whether package 'mrgsolve' can be installed ... OK
* used C++ compiler: 'g++.exe (GCC) 13.2.0'
* checking installed package size ... OK
* checking package directory ... OK
* checking for future file timestamps ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking code files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking whether startup messages can be suppressed ... OK
* checking use of S3 registration ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... [21s] OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking contents of 'data' directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking use of PKG_*FLAGS in Makefiles ... OK
* checking use of SHLIB_OPENMP_*FLAGS in Makefiles ... OK
* checking pragmas in C/C++ headers and code ... OK
* checking compilation flags used ... OK
* checking compiled code ... OK
* checking examples ... [16s] OK
* checking for unstated dependencies in 'tests' ... OK
* checking tests ... [21s] OK
  Running 'testthat.R' [21s]
* checking PDF version of manual ... [18s] OK
* checking HTML version of manual ... [20s] OK
* DONE
Status: 1 NOTE

@kylebaron kylebaron merged commit 683bf7b into main Sep 25, 2024
7 checks passed
@kylebaron kylebaron deleted the event-block branch September 25, 2024 12:31
@kylebaron kylebaron mentioned this pull request Oct 17, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants