Segfaultish issue when using ii=1 and ss=1 for ADVAN 2

[dpastoor]

Sorry for brevity on the road so no (computer) Internet will update with code examples when I'm back Monday wanted to put this as a placeholder.

Completely kills the R and rstudio processes

[kylebaron]

infusion?

[dpastoor]

I'll call it an 'accidental' infusion, example as shown below

The goal is the infusion should always be given over the course of an hour, so given an ii of 1, would essentially be an infusion.

library(mrgsolve)
library(dplyr)
code <- '
$PARAM CL=0.09, VC=1.1
$CMT GUT CENT
$ADVAN2
$MAIN
pred_CL = CL;
pred_V = VC;
'
model <- mread(code=code, model="onecmpt")

conc_range <- function(model, .amt, .ii) {
  e <- ev(amt=.amt,rate=.amt,ii=.ii,cmt=2,ss=1)
  result <- model %>% ev(e) %>%
    Req(CENT) %>% mrgsim(end=24,delta=0.1,digits=5) %>%
    filter(time > 0.5) %>%
    summarize(max_conc = max(CENT),
              min_conc = min(CENT)) %>% mutate(amt = .amt, ii = .ii) 
  return(result)
}

# will cause crash
conc_time(model, 100, 1) # found by setting default param of .ii=1 while testing function

# fine
conc_time(model, 100, 12)

Running the code interactively (outside function context) causes the same crash

defensive coding solution for now

conc_range <- function(model, .amt, .ii) {
 # ii should not be equal or less than infusion interval
    if (.ii <=1) {
        stop("ii cannot be less than or equal to the infusion interval")
    }
    ...
}

[kylebaron]

Thx; I'll look into it this weekend.

[kylebaron]

I don't think it's ADVAN4 per se; you should be able to do infusions that last longer than the dosing interval (with analytical equations or odes), but it looks like mrgsolve can't handle the ss flag when infusion duration is not strictly > ii. When I do the same sort of thing with ODE model, it just keeps working on something.

Bill Gillespie has ss equations for ADVAN2 and 4 that I'm not taking advantage of. I'll try to get those working but will need to look again at the ss algorithm for ODEs when ii > infusion duration. At least I'll enforce dosing interval greater than infusion interval when ss flag as limitation for now to prevent crash.

[dpastoor]

Sounds reasonable to me. FYI, I also did a benchmark of the closed form and ODE solutions and the timings were virtually identical (one difference was the ODE solution timings were much more variable, but mean and median times were within 0.5 ms of one another) - any reason for this? All the benchmarks I've done using python/julia/desolve+R see ~10x speedup of closed form vs the ODE implementations. Any idea why?

Could this be because it due to bottleneck of data exchange or some disk I/O?

I was using the test13.R comparison of ODE to analytical solutions for ADVAN2 examples for the benchmarks.

[kylebaron]

I saw the same thing: http://metrumresearchgroup.github.io/mrgsolve/advan4.html When the problem is small you don't see that much of a difference between analytical and ode. But you do see ~7x speedup when there are a lot of individuals simulated for a long time.

In the unit tests, I limited hmax to 0.1 and decreased the tolerances to 1E-12 in order to get the same answers between analytical and ode. That was slowing things down but obviously can't explain the 10x speedup you were seeing.

There is no disk IO for either ODE or analytical or data exchange differences: ADVAN2/4 are identical to the ODE under the hood, but for the analytical methods, it calls C++ algebraic code instead of the solver to advance the system from tfrom to tto. All of the other overhead is the same. The algebraic code was taken from Bill's bugs model library https://bitbucket.org/metrumrg/bugsmodellibrary/wiki/Home that Charles had translated into C++. It is a general poly exponential function that calculates amounts for all compartments handling all of the infusions etc. Maybe it could be speeded up if we're only doing bolus dosing and observing a single compartment?? Not sure ....

[kylebaron]

Just reporting back / thinking out loud after looking into this:

I don't totally have it nailed down, but I don't the problem is ss=1 and dur > ii in general. We can handle that and going back and checking we are properly tracking the number of infusions that are still active when ss has been satisfied.

It seems like the issue comes up when the infusion duration is some multiple of the dosing interval.

So this chokes:

ev(amt=100,rate=amt, ii=1)

because ii==dur

I was earlier looking at this:

ev(amt=100,rate=amt, ii=0.5)

worried that either (1) super-frequent dosing or (2) ii < dur , both with ss=1 was the causing the break.

But this is fine:

ev(amt=100,rate=amt, ii=0.51)

Hopefully this is just a quick book keeping fix and we won't need to enforce any restrictions on ii and dur here.

Regarding the segfaultish behavior with $ADVAN2/4: I'll also go in and figure out where that was going haywire and try to get that to exit more gracefully. Hopefully once the ss gets straightened out it won't cause the segfaulty response but better know why that is happening in case it ever gets that type of input again.

[kylebaron]

It turns out the problem was none of the above.

• SS algorithm for infusions was working fine
• The polyExp function to do the calculations for ADVAN2/4 was fine
• Things were crashing because the simulation loop started recursively calling the SS algorithm

The problem was this: when the infusion duration > dosing interval, there will be some infusions running at steady state that need to be turned off at the right time. We were taking this into account, but the calculation wasn't quite right when duration == interval.

But thank you @dpastoor for raising this issue: there was another problem in there that I discovered and resolved when straightening this out ... that code to schedule those infusion off events assumed that the steady state infusion dosing record was given at time==0. Ugh. Anyway, that is fixed now too.