Skip to content

Macros to determine the last record for an individual or a simulation #1327

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.

Sign up for GitHub

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

Merged
kylebaron merged 9 commits into from
Feb 25, 2026
Merged

Macros to determine the last record for an individual or a simulation #1327

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
6bab7f6
add macros for FINAL rows
kylebaron Feb 13, 2026
e25c711
make macros reserved
kylebaron Feb 13, 2026
3c01d76
FINAL_ROW macro tests; no modeled events after last time
kylebaron Feb 13, 2026
d4c5f2a
increment row counter past number of data rows
kylebaron Feb 13, 2026
693ae4f
refactor to continue for any record beyond maxtime
kylebaron Feb 13, 2026
a6e7561
finish tests
kylebaron Feb 13, 2026
4d765d3
refactor NNI
kylebaron Feb 14, 2026
3f2c68c
update gh issue/pr number in test header
kylebaron Feb 17, 2026
3068c9e
revert
kylebaron Feb 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion R/Aaaa.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,8 @@ Reserved_cvar <- c("SOLVERTIME","table","ETA","EPS", "AMT", "CMT",
"NEWIND", "DONE", "CFONSTOP", "DXDTZERO",
"CFONSTOP","INITSOLV","_F", "_R","_ALAG",
"SETINIT", "report", "_VARS_", "VARS",
"SS_ADVANCE", "END_OF_INFUSION", "CHECK_MODELED_INFUSIONS")
"SS_ADVANCE", "END_OF_INFUSION", "CHECK_MODELED_INFUSIONS",
"FINAL_ROW", "FINAL_IROW")

Reserved <- c("ID", "amt", "cmt", "ii", "ss", "evid",
"addl", "rate","time", Reserved_cvar,
Expand Down
4 changes: 4 additions & 0 deletions inst/base/modelheader.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,10 @@ typedef double capture;
#define THETA(a) THETA##a
// Should modeled infusion parameters get checked
#define CHECK_MODELED_INFUSIONS _check_modeled_infusions
// Check for final record of the simulation
#define FINAL_ROW (self.rown == (self.nrow-1))
// Check for final record of current individual
#define FINAL_IROW (self.irown == (self.inrow-1))

// NMVARS
#ifdef _MRGSOLVE_USING_NM_VARS_
Expand Down
60 changes: 57 additions & 3 deletions inst/maintenance/unit/test-rown-nrow.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,19 +234,33 @@ test_that("individual row counters work with PRED model", {

# https://github.com/metrumresearchgroup/mrgsolve/pull/1323
code_counter_update_on_output <- '
$preamble capture total1 = 0; capture total2 = 0;
$main self.mtime(12);
$preamble
capture total1 = 0; capture total2 = 0;
capture total3 = 0; capture total4 = 0;
capture total5 = 0;

$main
self.mtime(1.2); self.mtime(2.9);
self.mtime(11); self.mtime(12); self.mtime(24);

$table
if(self.rown+1==self.nrow) ++total1;
if(self.nid==self.idn+1 && self.irown+1 == self.inrow) ++total2;

if(FINAL_ROW) ++total3;
if(FINAL_IROW) ++total4;

if(self.irown+1==self.inrow) ++total5;

capture rown = self.rown;
capture nrow = self.nrow;
capture irown = self.irown;
capture inrow = self.inrow;
'

mod <- mcode("counter-update-on-output", code_counter_update_on_output)

test_that("row counters are only updated on output records gh-1323", {
mod <- mcode("counter-update-on-output", code_counter_update_on_output)
data <- expand.ev(amt = 0, cmt = 0, ID = 1:2)
out <- mrgsim(mod, data = data, end = 24, delta = 24)
expect_equal(nrow(out), 6)
Expand All @@ -263,7 +277,47 @@ test_that("row counters are only updated on output records gh-1323", {
expect_equal(out$total1[6], 1)
expect_equal(out$total2[6], 1);

# Macros
expect_equal(out$total1, out$total3)
expect_equal(out$total4, out$total4)
})

test_that("more row counter macros gh-1327", {
data1 <- ev(amt = 0, cmt = 1, ID = 1)
data2 <- ev(amt = 0, cmt = 1, ID = 2)

data1 <- expand_observations(data1, times = 1:5)
data2 <- expand_observations(data2, times = 1:8)

data <- rbind(data1, data2)
data$cmt <- 0

out <- mrgsim_df(mod, data = data)
expect_equal(nrow(out), nrow(data))

out1 <- subset(out, ID==1)

expect_equal(nrow(out1), 6)
expect_true(all(out1$inrow==6))
expect_equal(out1$rown, out1$irown)
expect_equal(out1$total4[6], 1)
expect_equal(out1$total4, out1$total5)
expect_equal(out1$irown, seq(6)-1)

expect_true(all(out1$total1==0))
expect_true(all(out1$total3==0))

out2 <- subset(out, ID==2)
expect_equal(nrow(out2), 9)
expect_true(all(out2$inrow==9))
expect_equal(out2$rown-6, out2$irown)
expect_equal(out2$total3[9], 1)
expect_equal(out2$total4, out2$total5)
expect_equal(out2$irown, seq(9)-1)

expect_true(all(out2$total4[1:8]==1))
})

rm(mod)
rm(code_test_rown_nrow, code_test_rown_nrow_pred)
rm(code_counter_update_on_output)
28 changes: 16 additions & 12 deletions src/devtran.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -427,7 +427,9 @@ Rcpp::List DEVTRAN(const Rcpp::List parin,

double tfrom = a[i].front()->time();
double tto = tfrom;
double maxtime = a[i].back()->time();

const double maxtime = a[i].back()->time();
const int NNI = dat.inrow(i);
Copy link
Collaborator Author

@kylebaron kylebaron Feb 17, 2026

Choose a reason for hiding this comment

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

The total row count for "this" individual.


for(int k=0; k < neta; ++k) prob.eta(k,eta(i,k));
for(int k=0; k < neps; ++k) prob.eps(k,eps(crow,k));
Expand Down Expand Up @@ -459,18 +461,22 @@ Rcpp::List DEVTRAN(const Rcpp::List parin,
Rcpp::checkUserInterrupt();
ic = prob.interrupt;
}

rec_ptr this_rec = a[i][j];

if(crow == NN) continue;
this_rec->id(id);
tto = this_rec->time();

// TODO: simplify
if(icrow==NNI || crow==NN || tto > maxtime) {
Copy link
Contributor

@kyleam kyleam Feb 17, 2026

Choose a reason for hiding this comment

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

Is it correct that this new tto > maxtime condition isn't an essential change, and it just avoids the downstream ->schedule call being a no-op? I removed the condition and the tests still pass.

Copy link
Collaborator Author

@kylebaron kylebaron Feb 17, 2026

Choose a reason for hiding this comment

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

Yeah; it's probably on the conservative side. I discovered that it was possible for some in consequential events happening after we hit the last official output record for a given individual. For example, we only calculate the end time of an infusion at the time when it starts, not when the dose is scheduled to happen (potentially in the future through addl) and I realized that we are allowing infusions to end after the last output record and the system is advancing to that time. I'm not totally sure what the right thing to do there. I'm guessing we could just decline to schedule the "infusion off" event if it would happen after maxtime. But I didn't want to get into that in this PR.

But you're right ... if maxtime (specific for an individual) is tied to the time when icrow==NNI, then we probably don't event need crow==NN, right? the NN test and the NNI test should both be true on the last record of the data set. I originally had this continue clause in there because crow references a row in the output matrix and you'd get bad behavior in case that count was ever off, so it's sort of fear-based decision :).

Could you test with just icrow==NNI?

      if(icrow==NNI || crow==NN || tto > maxtime) {
        continue;
      }

Copy link
Collaborator Author

@kylebaron kylebaron Feb 17, 2026

Choose a reason for hiding this comment

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

Also ... I did open an issue to at least review cases where we can be more careful of turning stuff off.
#1328

Copy link
Contributor

@kyleam kyleam Feb 17, 2026

Choose a reason for hiding this comment

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

Thanks for the details/explanation. Leaning conservative/safe is of course fine by me.

Could you test with just icrow==NNI?

Yes, the tests still all pass on my end with the condition reduced to if(icrow==NNI).

Copy link
Collaborator Author

@kylebaron kylebaron Feb 17, 2026

Choose a reason for hiding this comment

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

Thanks for checking that. I was almost going to push a commit that removed tto > maxtime, but think I'm going to leave it for now. I still need to fix the issue in #1324 and I think once that is behaving, I can write some check for correctness of the counts. At that point, will go back and see if we can simplify these checks, which I do think we should be able to do.

continue;
}

rec_ptr this_rec = a[i][j];

// Only update row counters on output records
if(this_rec->output()) {
prob.irown(icrow);
prob.rown(crow);
}

this_rec->id(id);
Copy link
Collaborator Author

@kylebaron kylebaron Feb 13, 2026

Choose a reason for hiding this comment

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

Moved this up for clarity.

}

if(prob.systemoff()) {
// This starts a loop that will finish the remaining records
Expand All @@ -483,8 +489,8 @@ Rcpp::List DEVTRAN(const Rcpp::List parin,
if(status==999) CRUMP("999 sent from the model.");
if(this_rec->output()) {
if(status==1) {
ans(crow,0) = this_rec->id();
ans(crow,1) = this_rec->time();
ans(crow,0) = id;
ans(crow,1) = tto;
for(unsigned int k=0; k < n_capture; ++k) {
Copy link
Collaborator Author

@kylebaron kylebaron Feb 17, 2026

Choose a reason for hiding this comment

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

Make consistent with main output writing code.

ans(crow,(k+capture_start)) = prob.capture(capture[k]);
}
Expand All @@ -511,9 +517,7 @@ Rcpp::List DEVTRAN(const Rcpp::List parin,
&prob
);
}

tto = this_rec->time();


double dt = (tto-tfrom)/(tfrom == 0.0 ? 1.0 : tfrom);

if((dt > 0.0) && (dt < mindt)) {
Expand Down
Loading