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

0 converted to small double/float when used in a power function #775

Merged
merged 8 commits into from
Aug 29, 2024

Conversation

mattfidler
Copy link
Member

Hello,

I found that 0 is converted to a small double (1.490116e-08) when used in the base of a power function inside an rxode2 model. This has caused large effects in my model. Particularly, when the amount of drug is 0 in the central compartment, I still see the "effect" compartment change, even though the derivative should be exactly 0 when the central compartment concentration is exactly 0. However, when I run the model with no drug added (i.e. central compartment is always exactly 0) the drug effect compartment changes 10-15% over the simulated 30 hours.

Here is a reproducible example where I've set some variables inside the rxode2 model definition. All of the variables that are 0^x should be 0 in the model results, however you can see in the output below that they are not.

  • 0^{1.2,0.67,1} are not 0, instead they are approximately (1.490116e-08) to the same values
  • 0 + 0 is giving the intended result: 0+0 = 0
  • 0/0 should be NaN but is 0, so it seems like the numerator is being treated as exactly 0 but the denominator is not (maybe a good thing in this context to avoid NAs)
  • x^0 = 0, which is expected so 0 is being treated as exactly 0 when it's the power value

Thanks!
Shayna

test_rxode_0 <- rxode2({
  fn(0) <- 0
  is_0_0 <- 0
  check_0_exp_1.2 <- 0^1.2
  check_near0_exp_1.2 <- (1.490116e-08)^1.2
  check_0_exp_0.67 <- 0^0.67
  check_near0_exp_0.67 <- (1.490116e-08)^0.67
  check_0_exp_1 <- 0^1
  check_near0_exp_near0 <- (1.490116e-08)^(1.490116e-08)
  check_near0_exp_0 <- (1.490116e-08)^0
  check_0_exp_0 <- 0^0
  check_0_exp_near0 <- 0^(1.490116e-08)
  check_0_div_0 <- 0/0
  check_0_plus_0 <- 0 + 0
  d/dt(fn) <- 1
})

test_events <-  et(seq(0,10, by = 1)) 

test_sim <- rxSolve(test_rxode_0, events = test_events)

test_sim %>% data.frame()

    time check_0_exp_1.2 check_near0_exp_1.2 check_0_exp_0.67 check_near0_exp_0.67 check_0_exp_1
1    0    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
2    1    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
3    2    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
4    3    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
5    4    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
6    5    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08
  check_near0_exp_near0 check_near0_exp_0 check_0_exp_0 check_0_exp_near0 check_0_div_0 check_0_plus_0 fn
1             0.9999997                 1             1         0.9999997             0              0  0
2             0.9999997                 1             1         0.9999997             0              0  1
3             0.9999997                 1             1         0.9999997             0              0  2
4             0.9999997                 1             1         0.9999997             0              0  3
5             0.9999997                 1             1         0.9999997             0              0  4
6             0.9999997                 1             1         0.9999997             0              0  5

@mattfidler
Copy link
Member

Hi @sstein93

You are right. This is about avoiding NaN issues.

You can turn this off by rxSolve(..., safeZero=FALSE)

@mattfidler
Copy link
Member

I have thought about refactoring this code for a while. In theory, different models should be generated based on if you want this protection or not (since it adds a bit of overhead). Then the powers and the like could probably be only fixed if there is an issue (though that would take a bit more overhead...)

@sstein93
Copy link
Author

Thank you @mattfidler!

@sstein93
Copy link
Author

@mattfidler I just tried the example above and adding safeZero = FALSE didn't seem to have any effect on the output. I also added a check for log(0) and dividing by 0.

Do you know if this is due to a specific version of an R or C dependency?

Thanks,
Shayna

test_rxode_0 <- rxode2({
  fn(0) <- 0
  is_0_0 <- 0
  check_0_exp_1.2 <- 0^1.2
  check_near0_exp_1.2 <- (1.490116e-08)^1.2
  check_0_exp_0.67 <- 0^0.67
  check_near0_exp_0.67 <- (1.490116e-08)^0.67
  check_0_exp_1 <- 0^1
  check_near0_exp_near0 <- (1.490116e-08)^(1.490116e-08)
  check_near0_exp_0 <- (1.490116e-08)^0
  check_0_exp_0 <- 0^0
  check_0_exp_near0 <- 0^(1.490116e-08)
  check_0_div_0 <- 0/0
  check_0_plus_0 <- 0 + 0
  check_log_0 <- log(0)
  check_divide_0 <- 1/0
  d/dt(fn) <- 1
})

test_events <-  et(seq(0,5, by = 1)) 

test_sim <- rxSolve(test_rxode_0, events = test_events,safeZero = FALSE)

  time check_0_exp_1.2 check_near0_exp_1.2 check_0_exp_0.67 check_near0_exp_0.67 check_0_exp_1 check_near0_exp_near0
1    0    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
2    1    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
3    2    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
4    3    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
5    4    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
6    5    4.053817e-10        4.053817e-10     5.702397e-06         5.702397e-06  1.490116e-08             0.9999997
  check_near0_exp_0 check_0_exp_0 check_0_exp_near0 check_0_div_0 check_0_plus_0 check_log_0 check_divide_0 fn
1                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  0
2                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  1
3                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  2
4                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  3
5                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  4
6                 1             1         0.9999997             0              0   -36.04365     4.5036e+15  5

Here is my sessionInfo in case that is helpful:

> sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.6.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rxode2_2.1.3.9000

loaded via a namespace (and not attached):
 [1] gtable_0.3.5         crayon_1.5.3         dplyr_1.1.4          compiler_4.3.2       tidyselect_1.2.1    
 [6] Rcpp_1.0.13          scales_1.3.0         fastmap_1.2.0        lattice_0.22-6       ggplot2_3.5.1       
[11] R6_2.5.1             generics_0.1.3       knitr_1.48           backports_1.5.0      checkmate_2.3.1     
[16] tibble_3.2.1         munsell_0.5.1        stringfish_0.16.0    pillar_1.9.0         rlang_1.1.4         
[21] utf8_1.2.4           cachem_1.1.0         xfun_0.46            sys_3.4.2            RcppParallel_5.1.7  
[26] memoise_2.0.1        cli_3.6.3            magrittr_2.0.3       digest_0.6.36        grid_4.3.2          
[31] rstudioapi_0.15.0    rxode2ll_2.0.11.9000 lotri_0.5.0.9000     lifecycle_1.0.4      nlme_3.1-164        
[36] dparser_1.3.1-12     vctrs_0.6.5          glue_1.7.0           data.table_1.15.4    PreciseSums_0.7     
[41] fansi_1.0.6          RApiSerialize_0.1.3  colorspace_2.1-0     tools_4.3.2          qs_0.26.3           
[46] pkgconfig_2.0.3     

@mattfidler
Copy link
Member

This is likely based on the version of rxode2 you are using. I plan
to refactor this to help your issue now.

@mattfidler
Copy link
Member

Hi @sstein93,

The pull request this is being changed into probably behaves a bit
better.

library(rxode2)
#> rxode2 2.1.3.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
test_rxode_0 <- rxode2({
  fn(0) <- 0
  is_0_0 <- 0
  check_0_exp_1.2 <- 0^1.2
  check_near0_exp_1.2 <- (1.490116e-08)^1.2
  check_0_exp_0.67 <- 0^0.67
  check_near0_exp_0.67 <- (1.490116e-08)^0.67
  check_0_exp_1 <- 0^1
  check_near0_exp_near0 <- (1.490116e-08)^(1.490116e-08)
  check_near0_exp_0 <- (1.490116e-08)^0
  check_0_exp_0 <- 0^0
  check_0_exp_near0 <- 0^(1.490116e-08)
  check_0_div_0 <- 0/0
  check_0_plus_0 <- 0 + 0
  d/dt(fn) <- 1
})

test_events <-  et(seq(0,10, by = 1))

test_sim <- rxSolve(test_rxode_0, events = test_events, returnType="data.frame",
        safeZero=TRUE, safePow=TRUE, safeLog=TRUE, useStdPow=FALSE)

print(test_sim)
#>    time check_0_exp_1.2 check_near0_exp_1.2 check_0_exp_0.67
#> 1     0               0        4.053817e-10                0
#> 2     1               0        4.053817e-10                0
#> 3     2               0        4.053817e-10                0
#> 4     3               0        4.053817e-10                0
#> 5     4               0        4.053817e-10                0
#> 6     5               0        4.053817e-10                0
#> 7     6               0        4.053817e-10                0
#> 8     7               0        4.053817e-10                0
#> 9     8               0        4.053817e-10                0
#> 10    9               0        4.053817e-10                0
#> 11   10               0        4.053817e-10                0
#>    check_near0_exp_0.67 check_0_exp_1 check_near0_exp_near0 check_near0_exp_0
#> 1          5.702397e-06             0             0.9999997                 1
#> 2          5.702397e-06             0             0.9999997                 1
#> 3          5.702397e-06             0             0.9999997                 1
#> 4          5.702397e-06             0             0.9999997                 1
#> 5          5.702397e-06             0             0.9999997                 1
#> 6          5.702397e-06             0             0.9999997                 1
#> 7          5.702397e-06             0             0.9999997                 1
#> 8          5.702397e-06             0             0.9999997                 1
#> 9          5.702397e-06             0             0.9999997                 1
#> 10         5.702397e-06             0             0.9999997                 1
#> 11         5.702397e-06             0             0.9999997                 1
#>    check_0_exp_0 check_0_exp_near0 check_0_div_0 check_0_plus_0 fn
#> 1              1                 0             0              0  0
#> 2              1                 0             0              0  1
#> 3              1                 0             0              0  2
#> 4              1                 0             0              0  3
#> 5              1                 0             0              0  4
#> 6              1                 0             0              0  5
#> 7              1                 0             0              0  6
#> 8              1                 0             0              0  7
#> 9              1                 0             0              0  8
#> 10             1                 0             0              0  9
#> 11             1                 0             0              0 10

test_sim <- rxSolve(test_rxode_0, events = test_events, returnType="data.frame",
    safeZero=FALSE, safePow=FALSE, safeLog=FALSE, useStdPow=TRUE)

# This tests some of the differences between the solves:

m <- rxode2({
  nan0 = 0/0
  nanLog = log(-1)
  nanPow = 0^-1
  rPow = 0^-1.5
})

et <- et(1)

# No safe zero in place, use C's pow() function
rxSolve(m, et, safeZero=FALSE, safePow=FALSE, safeLog=FALSE, useStdPow = TRUE, returnType="data.frame")
#>   time nan0 nanLog nanPow rPow
#> 1    1  NaN    NaN    Inf  Inf

# No safe zero in place, use R's pow() function
rxSolve(m, et, safeZero=FALSE, safePow=FALSE, safeLog=FALSE, useStdPow = FALSE, returnType="data.frame")
#>   time nan0 nanLog nanPow rPow
#> 1    1  NaN    NaN    Inf  Inf

# Safe zero for everything
rxSolve(m, et, safeZero=TRUE, safePow=TRUE, safeLog=TRUE, useStdPow = FALSE, returnType="data.frame")
#>   time nan0    nanLog     nanPow         rPow
#> 1    1    0 -36.04365 4.5036e+15 3.022315e+23

Created on 2024-08-28 with reprex v2.1.1

The rules are updated to:

  • when safeZero=TRUE and the denominator of a division expression is
    zero, use the Machine's small number/eps (you can see this value
    with .Machine$double.eps)

  • when saveLog=TRUE and the x in the log(x) is less than or equal
    to zero, change this to log(eps)

  • when safePow=TRUE and the expression x^y has a zero for x and
    a negative number for y replace x with eps.

  • when useStdPow=TRUE use C's internal pow() function. When
    useStdPow=FALSE use R's internal R_pow() function. You can see
    the difference in the two implementations when taking 0^-0.5.

If this successfully passes the checks, I will merge into the main
branch and then you should be able to get this from the development
R-universe repository.

I also plan to release at the end of September 2024 due to a CRAN
request.

@mattfidler
Copy link
Member

Thanks for reporting @sstein93

@mattfidler mattfidler merged commit b1bf8b0 into main Aug 29, 2024
5 of 7 checks passed
@mattfidler mattfidler deleted the 775-solveZero branch August 29, 2024 01:29
@sstein93
Copy link
Author

sstein93 commented Sep 4, 2024

Thanks @mattfidler !!

@mattfidler
Copy link
Member

No problem.

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