Correctly handle EVID4 doses at steady state with lag time

inst/maintenance/unit/test-alag.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2013 - 2019  Metrum Research Group
+# Copyright (C) 2013 - 2025  Metrum Research Group
 #
 # This file is part of mrgsolve.
 #
@@ -97,3 +97,21 @@ test_that("ss dose with lag time issue-484", {
   out2 <- mrgsim_d(mod,data2) %>% slice(-c(1,2))
   expect_identical(out1,out2)
 })
+
+test_that("EVID 1 or 4 at SS with lag time are identical", {
+  ev1 <- ev(amt = 100, ii = 24, ss = 1, LAG = 10)
+  ev4 <- mutate(ev1, evid = 4)
+
+  out1 <- mrgsim_e(mod, ev1, recsort = 3, output = "df")
+  out4 <- mrgsim_e(mod, ev4, recsort = 3, output = "df")
+
+  expect_identical(out1, out4)
+
+  ev1a <- mutate(ev1, addl = 3, LAG = 5)
+  ev4a <- mutate(ev4, addl = 3, LAG = 5)
+
+  out1a <- mrgsim_e(mod, ev1a, recsort = 3, output = "df")
+  out4a <- mrgsim_e(mod, ev4a, recsort = 3, output = "df")
+
+  expect_identical(out1a, out4a)
+})

mrgsolve.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 0063b92e-a050-4fd7-986e-be62f57eb1b1
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

src/datarecord.cpp

@@ -191,7 +191,8 @@ void datarecord::implement(odeproblem* prob) {
   case 8: // replace
     prob->y(eq_n, Amt);
     break;
-  case 4:
+  case 4: // dose and reset; reset only if non-ss
+    if(!Armed) break;
     if(this->ss()==0) {
       for(int i=0; i < prob->neq(); ++i) {
         prob->y(i,0.0);
@@ -200,7 +201,6 @@ void datarecord::implement(odeproblem* prob) {
       }
       prob->init_call(Time);
     }
-    if(!Armed) break;
     if(Rate > 0) {
       this->evid(5);
     } else {

src/devtran.cpp

@@ -548,19 +548,23 @@ Rcpp::List DEVTRAN(const Rcpp::List parin,
         // Checking 
         if(!this_rec->is_lagged()) {
 
-          if(prob.alag(this_cmtn) > mindt && this_rec->is_dose()) { // there is a valid lagtime
-            
+          // there is a valid lag time
+          if(prob.alag(this_cmtn) > mindt && this_rec->is_dose()) {
             if(this_rec->ss() > 0) {
               this_rec->steady(&prob, a[i], solver);
               tfrom = tto;
-              this_rec->ss(0);
             }
+            // We already advanced to ss
+            // Lagged dose and all subsequent should be vanilla EVID=1 doses
+            this_rec->ss(0);
             rec_ptr newev = NEWREC(*this_rec);
+            if(newev->evid()==4) {
+              newev->evid(1);
+            }
             newev->pos(__ALAG_POS);
             newev->phantom_rec();
             newev->lagged();
             newev->time(this_rec->time() + prob.alag(this_cmtn));
-            newev->ss(0);
             insert_record(a[i], j, newev, put_ev_first);
             newev->schedule(a[i], maxtime, put_ev_first, NN, prob.alag(this_cmtn));
             this_rec->unarm();