search

nlmixrdevelopment / RxODE

RxODE is an R package that facilitates easy simulations in R
https://nlmixrdevelopment.github.io/RxODE/
GNU General Public License v3.0
54 stars 14 forks source link

Error in rxSolveSEXP(...) could not allocate enough memory for solving #459

Closed levenc closed 3 years ago

levenc commented 3 years ago

Hello, The following error appeared when applying rxSolve() to a dataset of 19,000 rows (10,000 individuals, with 19 rows per patient dataset):

Error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits,  : 
  could not allocate enough memory for solving

My system: Dell XPS 13 9360, 16 GB RAM, 7th generation Intel Core i7 processor.

Surprisingly, whatever the memory available on the system, the error always occurs for the same number of individuals: no problem for n = 7,600 individuals, error for 7,700 individuals. I tested the same code today on a server with twice as much RAM, and the values that trigger the error were identical (with the same version of R, same version of RxODE, same GNU/linux distribution).

Here is a code snippet that reproduces the error:

# libraries
library(lotri)
library(RxODE)
library(data.table)
library(mvtnorm)

# number of samples
n_samples <- 7.7e3

# RxODE model
ppk_model   = RxODE({
  centr(0)   = 0;
  # Time lag from dose infuser to patient
  TVLAGTIME  = THETA_LAGTIME
  LAGTIME    = TVLAGTIME*exp(ETA_LAGTIME)
  # -------------------------------------
  TVCl  = THETA_Cl*(TBW/70);
  TVV1  = THETA_V1;
  TVV2  = THETA_V2;
  TVQ   = THETA_Q;
  Cl    = TVCl*exp(ETA_Cl + KAPPA_Cl);
  Cli   = TVCl*exp(ETA_Cl);
  V1    = TVV1*exp(ETA_V1 + KAPPA_V1);
  V1i   = TVV1*exp(ETA_V1);
  V2    = TVV2*exp(ETA_V2);
  Q     = TVQ *exp(ETA_Q);
  ke    = Cl/V1;
  k12   = Q/V1;
  k21   = Q/V2;
  Cc    = centr/V1;
  d/dt(centr)  = - ke*centr - k12*centr + k21*periph;
  d/dt(periph) =            + k12*centr - k21*periph;
  lag(centr)   =   LAGTIME;
  d/dt(AUC)    =   Cc;
})

# patient
datapat       <- data.frame(ID=6,TIME=c(0.0,8,16,24,32,40,47.9,48,
                            56,64,71.9,72,80,88,95.9,96,104,112,120),
                            DV=c(NA,NA,NA,NA,NA,NA,63.0,NA,NA,NA,
                            19.7,NA,NA,NA,31,NA,NA,NA,16),
                            AMT=c(4000,4000,4000,4000,4000,4000,NA,
                            4000,4000,4000,NA,4000,4000,4000,NA,4000,
                            4000,4000,NA),
                            DUR=c(8,8,8,8,8,8,NA,8,8,8,NA,8,8,8,NA,8,
                            8,8,NA),
                            EVID=c(1,1,1,1,1,1,0,1,1,1,0,1,1,1,0,1,1,
                            1,0),OCC=c(1,1,1,1,1,1,1,2,2,2,2,3,
                            3,3,3,4,4,4,4,TBW=92,))
# theta vector
theta = rbind(c(THETA_LAGTIME=0.07,THETA_Cl=17.1,THETA_V1=7.2,THETA_V2=17.8,
          THETA_Q=52.0))

# omega matrix
omega = lotri({ETA_LAGTIME + ETA_Cl + ETA_V1 + ETA_V2 + ETA_Q + KAPPA_Cl_1 +
    KAPPA_V1_1 + KAPPA_Cl_2 + KAPPA_V1_2 + KAPPA_Cl_3 + KAPPA_V1_3 +
    KAPPA_Cl_4 + KAPPA_V1_4 ~
    c(0.1747673,
      0.00     ,  0.08507985,
      0.00     ,  0.00      ,  0.0673745,
      0.00     ,  0.00      ,  0.00     ,  0.429067,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.2247455,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.1934626,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.05783106,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.1934626,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.00    ,   0.05783106,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.00    ,   0.00    ,    0.1934626,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.05783106,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.1934626,
      0.00     ,  0.00      ,  0.00     ,  0.00    ,  0.00     ,  0.00    ,  0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.00    ,   0.05783106)})

# simulate random effects
eta_sim           <- rmvnorm(n_samples,mean=rep(0,ncol(omega)),
                                      sigma=omega)
colnames(eta_sim) <- attr(omega,"dimnames")[[1]]
eta_dt            <- data.table(eta_sim)
eta_dt[,ID:=1:n_samples]

dat_dt <- data.table(datapat)
dat_dt <- dat_dt[rep(dat_dt[,.I],n_samples)]         # n_samples patient records
dat_dt[,ID:=rep(1:n_samples,1,each=nrow(datapat))]

# bind random effects to patient data.table
dat_dt <- dat_dt[eta_dt,on = .(ID = ID), roll = TRUE]

dat_dt[OCC == 1, c("KAPPA_Cl", "KAPPA_V1") := list(KAPPA_Cl_1,KAPPA_V1_1), by = ID]
dat_dt[OCC == 2, c("KAPPA_Cl", "KAPPA_V1") := list(KAPPA_Cl_2,KAPPA_V1_2), by = ID]
dat_dt[OCC == 3, c("KAPPA_Cl", "KAPPA_V1") := list(KAPPA_Cl_3,KAPPA_V1_3), by = ID]
dat_dt[OCC == 4, c("KAPPA_Cl", "KAPPA_V1") := list(KAPPA_Cl_4,KAPPA_V1_4), by = ID]

dat_dt <- dat_dt[,!c("KAPPA_Cl_1","KAPPA_V1_1","KAPPA_Cl_2","KAPPA_V1_2",
                     "KAPPA_Cl_3","KAPPA_V1_3","KAPPA_Cl_4","KAPPA_V1_4")]

solved_model <- RxODE::rxSolve(ppk_model,cbind(theta,eta_dt),dat_dt)

─ Session info ──────────────────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 4.1.1 (2021-08-10)
 os       Ubuntu 20.04.3 LTS          
 system   x86_64, linux-gnu           
 ui       RStudio                     
 language fr_FR                       
 collate  fr_FR.UTF-8                 
 ctype    fr_FR.UTF-8                 
 tz       Europe/Paris                
 date     2021-09-21                  

─ Packages ──────────────────────────────────────────────────────────────────────────────────────────
 package       * version date       lib source        
 backports       1.2.1   2020-12-09 [1] CRAN (R 4.1.0)
 cachem          1.0.6   2021-08-19 [1] CRAN (R 4.1.1)
 callr           3.7.0   2021-04-20 [1] CRAN (R 4.1.0)
 checkmate       2.0.0   2020-02-06 [1] CRAN (R 4.1.0)
 cli             3.0.1   2021-07-17 [1] CRAN (R 4.1.0)
 colorspace      2.0-2   2021-06-24 [1] CRAN (R 4.1.0)
 crayon          1.4.1   2021-02-08 [1] CRAN (R 4.1.0)
 data.table    * 1.14.0  2021-02-21 [1] CRAN (R 4.1.0)
 desc            1.3.0   2021-03-05 [1] CRAN (R 4.1.0)
 devtools        2.4.2   2021-06-07 [1] CRAN (R 4.1.0)
 digest          0.6.27  2020-10-24 [1] CRAN (R 4.1.0)
 dparser         1.3.1-4 2021-04-07 [1] CRAN (R 4.1.0)
 dplyr           1.0.7   2021-06-18 [1] CRAN (R 4.1.0)
 ellipsis        0.3.2   2021-04-29 [1] CRAN (R 4.1.0)
 fansi           0.5.0   2021-05-25 [1] CRAN (R 4.1.0)
 fastmap         1.1.0   2021-01-25 [1] CRAN (R 4.1.0)
 fs              1.5.0   2020-07-31 [1] CRAN (R 4.1.0)
 generics        0.1.0   2020-10-31 [1] CRAN (R 4.1.0)
 ggplot2         3.3.5   2021-06-25 [1] CRAN (R 4.1.0)
 glue            1.4.2   2020-08-27 [1] CRAN (R 4.1.0)
 gtable          0.3.0   2019-03-25 [1] CRAN (R 4.1.0)
 knitr           1.34    2021-09-09 [1] CRAN (R 4.1.1)
 lifecycle       1.0.0   2021-02-15 [1] CRAN (R 4.1.0)
 lotri         * 0.3.1   2021-01-05 [1] CRAN (R 4.1.0)
 magrittr        2.0.1   2020-11-17 [1] CRAN (R 4.1.0)
 memoise         2.0.0   2021-01-26 [1] CRAN (R 4.1.0)
 munsell         0.5.0   2018-06-12 [1] CRAN (R 4.1.0)
 mvtnorm       * 1.1-2   2021-06-07 [1] CRAN (R 4.1.0)
 pillar          1.6.2   2021-07-29 [1] CRAN (R 4.1.0)
 pkgbuild        1.2.0   2020-12-15 [1] CRAN (R 4.1.0)
 pkgconfig       2.0.3   2019-09-22 [1] CRAN (R 4.1.0)
 pkgload         1.2.2   2021-09-11 [1] CRAN (R 4.1.1)
 PreciseSums     0.4     2020-07-10 [1] CRAN (R 4.1.0)
 prettyunits     1.1.1   2020-01-24 [1] CRAN (R 4.1.0)
 processx        3.5.2   2021-04-30 [1] CRAN (R 4.1.0)
 ps              1.6.0   2021-02-28 [1] CRAN (R 4.1.0)
 purrr           0.3.4   2020-04-17 [1] CRAN (R 4.1.0)
 qs              0.25.1  2021-07-25 [1] CRAN (R 4.1.0)
 R6              2.5.1   2021-08-19 [1] CRAN (R 4.1.1)
 RApiSerialize   0.1.0   2014-04-19 [1] CRAN (R 4.1.0)
 Rcpp            1.0.7   2021-07-07 [1] CRAN (R 4.1.0)
 RcppParallel    5.1.4   2021-05-04 [1] CRAN (R 4.1.0)
 remotes         2.4.0   2021-06-02 [1] CRAN (R 4.1.0)
 rlang           0.4.11  2021-04-30 [1] CRAN (R 4.1.0)
 rprojroot       2.0.2   2020-11-15 [1] CRAN (R 4.1.0)
 rstudioapi      0.13    2020-11-12 [1] CRAN (R 4.1.0)
 RxODE         * 1.1.1   2021-09-03 [1] CRAN (R 4.1.1)
 scales          1.1.1   2020-05-11 [1] CRAN (R 4.1.0)
 sessioninfo     1.1.1   2018-11-05 [1] CRAN (R 4.1.0)
 stringfish      0.15.2  2021-07-24 [1] CRAN (R 4.1.0)
 sys             3.4     2020-07-23 [1] CRAN (R 4.1.0)
 testthat        3.0.4   2021-07-01 [1] CRAN (R 4.1.0)
 tibble          3.1.4   2021-08-25 [1] CRAN (R 4.1.1)
 tidyselect      1.1.1   2021-04-30 [1] CRAN (R 4.1.0)
 usethis         2.0.1   2021-02-10 [1] CRAN (R 4.1.0)
 utf8            1.2.2   2021-07-24 [1] CRAN (R 4.1.0)
 vctrs           0.3.8   2021-04-29 [1] CRAN (R 4.1.0)
 withr           2.4.2   2021-04-18 [1] CRAN (R 4.1.0)
 xfun            0.26    2021-09-14 [1] CRAN (R 4.1.1)

On the other hand, updating the parameters of a model for 100,000 patients with the following syntax is not a problem :

solved_model$params  <- cbind(theta,eta_dt,row.names=NULL)

Unfortunately, I am not using this approach here because I use the patient data set to input inter-occasion variability as time-varying covariates (I use RxODE as a simulation engine for a dosage individualization R package I'm working on, and for now the randoms effects are handled separately in R).

Also, and unrelated to this error, I had not had the opportunity to thank you for issue 185: thank you for your help.

levenc commented 3 years ago

Some additional information after my tests this morning.

Using a single vector of parameters for all individuals, I get no errors. 

solved_model <- RxODE::rxSolve(ppk_model,cbind(theta,eta_dt[1,!c("ID")]),dat_dt)

However, when I try to include all the individual parameters, the same error occurs. 

solved_model$params <- cbind(theta,eta_dt)
Error in `$<-.rxSolve`(`*tmp*`, params, value = list(THETA_LAGTIME = c(0.07,  : 
  could not allocate enough memory for solving

One clarification (this may be a problem on my side), when I use a single vector of parameters with the model of the original post, the solving fails with this message:

solved_model <- RxODE::rxSolve(ppk_model,cbind(theta,eta_dt[1,!c("ID")]),dat_dt)
[====|====|====|====|====|====|====|====|====|====] 0:00:00 

Error : 'alag(.)'/'rate(.)'/'dur(.)' cannot depend on the state values

I had to comment out the lines that relate to LAGTIME in the model to get it to work.

mattfidler commented 3 years ago

Hi @levenc

I cannot solve a problem of running out of memory very easily. But I can give you a possible work-around. Simply solve the first group of subjects and then solve the second group of subjects. In between save the files and then put them together in the end.

The memory that RxODE currently requires is the:

Some memory could be spared by not duplicating the data, but this would also remove the ability to solve the system in parallel. In R, the data structures like data.frames are actually not thread safe. This would also require a major rewrite of the solving engine to support both duplicate data structures for solving the system in a threaded manner as well as another data structure that takes the data in-place. While I think it is valuable, it is not the top of my todo list currently.

I'm unsure if this is helpful, but it at least explains the current state of RxODE.

mattfidler commented 3 years ago

Note you could also use:

rxSolve(..., returnType="data.table")

Since you seem to be using data.table elsewhere. It will also save memory and slightly speed up the solving.

mattfidler commented 3 years ago

Also note, there is some rudimentary support for inter-occasion varaibility built in to the latest RxODE. It itsn't as well tested, but you can see the example here:

https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-nesting.html?q=occasion#nesting-variability

levenc commented 3 years ago

Thank you @mattfidler for this detailed answer, I understand the situation better.

To get around the problem, I will implement the option you suggest: divide the subject group into smaller sets, solve for each set, and then combine them afterwards.

Also thanks for the note on data.table, that's good to know. This is the first time I use data.table for one of my projects, but certainly not the last; the performance and syntax make me want to forget about data.frames.