R Under development (unstable) (2024-06-17 r86768 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(epigrowthfit)
> options(warn = 2L, error = if (interactive()) recover)
> 
> src <- "src"
> src. <- system.file("tests", src, package = "epigrowthfit", mustWork = TRUE)
> file.copy(c(src, src.), tempdir(), recursive = TRUE)
[1] TRUE TRUE
> setwd(file.path(tempdir(), src))
> Sys.setenv(R_TESTS = "") # startup.Rs does not exist here
> 
> dll <- "test"
> cpp <- paste0(dll, ".cpp")
> TMB::compile(cpp)
using C++ compiler: 'g++.exe (GCC) 13.2.0'
make[1]: Entering directory '/d/temp/RtmpCGpfGE/src'
g++ -std=gnu++17  -I"D:/RCompile/recent/R/include" -DNDEBUG -I"D:/RCompile/CRANpkg/lib/4.5/TMB/include" -I"D:/temp/Rtmp8uMMHd/RLIBS_ffac61196da9/RcppEigen/include"  -DTMB_SAFEBOUNDS -DTMB_EIGEN_DISABLE_WARNINGS -DLIB_UNLOAD=R_unload_test  -DTMB_LIB_INIT=R_init_test  -DCPPAD_FRAMEWORK    -I"d:/rtools44/x86_64-w64-mingw32.static.posix/include"     -pedantic -O2 -Wall  -mfpmath=sse -msse2 -mstackrealign  -c test.cpp -o test.o
g++ -std=gnu++17 -shared -s -static-libgcc -o test.dll tmp.def test.o -Ld:/rtools44/x86_64-w64-mingw32.static.posix/lib/x64 -Ld:/rtools44/x86_64-w64-mingw32.static.posix/lib -LD:/RCompile/recent/R/bin/x64 -lR
make[1]: Leaving directory '/d/temp/RtmpCGpfGE/src'
make[1]: Entering directory '/d/temp/RtmpCGpfGE/src'
make[1]: Leaving directory '/d/temp/RtmpCGpfGE/src'
[1] 0
> dyn.load(paste0(dll, .Platform[["dynlib.ext"]]))
> 
> getReport <- local({
+ 	tt <- readLines(cpp)
+ 	tt <- paste(tt, collapse = "\n")
+ 	tt <- chartr("\t\r\n", "   ", tt)
+ 	tt <- sub("^.*?enum +test *\\{(.*?)\\} *;.*$", "\\1", tt)
+ 	tt <- gsub(" ", "", tt)
+ 	enum. <- strsplit(tt, ",")[[1L]]
+ 
+ 	function(enum, ...)
+ 		TMB::MakeADFun(data = list(flag = match(enum, enum., 0L) - 1L, ...),
+ 		               parameters = list(),
+ 		               type = "Fun",
+ 		               checkParameterOrder = FALSE,
+ 		               DLL = dll)[["report"]]()[["ans"]]
+ })
> 
> mvlgamma <- function(x, p)
+     0.25 * p * (p - 1) * log(pi) + rowSums(lgamma(outer(x, seq.int(0, by = 0.5, length.out = p), `-`)))
> 
> dlkj <- function(x, eta, log = FALSE) {
+ 	n <- 0.5 * (1 + sqrt(1 + 8 * length(x)))
+ 	R <- diag(n)
+ 	R[upper.tri(R)] <- x
+ 	log.ans <- (eta - 1) * (-sum(log(colSums(R * R))))
+ 	if (log) log.ans else exp(log.ans)
+ }
> 
> dwishart <- function(x, df, scale, log = FALSE) {
+ 	n <- 0.5 * (-1 + sqrt(1 + 8 * length(x)))
+ 	X <- theta2cov(x)
+ 	S <- theta2cov(scale)
+ 	log.ans <- -0.5 * (df * log(det(S)) + (-df + n + 1) * log(det(X)) + n * df * log(2) + 2 * mvlgamma(0.5 * df, n) + sum(diag(solve(S, X))))
+ 	if (log) log.ans else exp(log.ans)
+ }
> 
> dinvwishart <- function(x, df, scale, log = FALSE) {
+ 	n <- 0.5 * (-1 + sqrt(1 + 8 * length(x)))
+ 	X <- theta2cov(x)
+ 	S <- theta2cov(scale)
+ 	log.ans <- -0.5 * (-df * log(det(S)) + (df + n + 1) * log(det(X)) + n * df * log(2) + 2 * mvlgamma(0.5 * df, n) + sum(diag(solve(X, S))))
+ 	if (log) log.ans else exp(log.ans)
+ }
> 
> 
> ## list_of_vectors_t
> x <- list(rnorm(10L), seq_len(5L), TRUE, double(0L))
> ans <- getReport("list_of_vectors_t", x = x)
> stopifnot(identical(ans, lapply(x, as.double)))
> 
> ## is_NA_real_
> x <- c(0, NA, NaN, Inf)
> ans <- getReport("is_NA_real_", x = x)
> stopifnot(identical(ans, c(0, 1, 0, 0)))
> 
> ## is_finite
> x <- c(0, NA, NaN, Inf)
> ans <- getReport("is_finite", x = x)
> stopifnot(identical(ans, c(1, 0, 0, 0)))
> 
> ## logspace_diff
> log.x <- log(cumsum(rlnorm(10L)))
> ans <- getReport("logspace_diff", log_x = log.x)
> stopifnot(all.equal(ans, log(diff(exp(log.x)))))
> 
> ## mvlgamma
> x <- seq.int(1, 12, by = 1)
> p1 <- 1L
> ans1 <- getReport("mvlgamma", x = x, p = p1)
> stopifnot(all.equal(ans1, mvlgamma(x, p1)))
> p2 <- 4L
> ans2 <- getReport("mvlgamma", x = x, p = p2)
> stopifnot(all.equal(ans2, mvlgamma(x, p2)))
> 
> ## dlkj
> n <- 4L
> x <- rnorm(0.5 * n * (n - 1L))
> eta <- 2
> ans <- getReport("dlkj", x = x, eta = eta, give_log = 1L)
> stopifnot(all.equal(ans, dlkj(x, eta, TRUE)))
> 
> ## d(inv)?wishart
> n <- 4L
> x <- rnorm(0.5 * n * (n + 1L))
> df <- 8
> scale <- rnorm(length(x))
> ans1 <- getReport(   "dwishart", x = x, df = df, scale = scale, give_log = 1L)
> stopifnot(all.equal(ans1,    dwishart(x, df, scale, TRUE)))
> ans2 <- getReport("dinvwishart", x = x, df = df, scale = scale, give_log = 1L)
> stopifnot(all.equal(ans2, dinvwishart(x, df, scale, TRUE)))
> 
> ## dpois_robust
> log.lambda <- seq.int(0, 10, by = 1)
> x <- rpois(log.lambda, lambda = exp(log.lambda))
> ans <- getReport("dpois_robust", x = x, log_lambda = log.lambda, give_log = TRUE)
> stopifnot(all.equal(ans, dpois(x, lambda = exp(log.lambda), log = TRUE)))
> 
> ## rnbinom_robust
> log.mu <- log(100)
> log.size <- log(50)
> n <- 1e+06L
> set.seed(10235L)
> ans <- getReport("rnbinom_robust", log_mu = log.mu, log_size = log.size, n = n)
> mm <- seq.int(min(ans), max(ans))
> freq <- as.integer(table(factor(ans, levels = mm)))
> dens <- dnbinom(mm, mu = exp(log.mu), size = exp(log.size))
> stopifnot(all.equal(freq / n, dens, tolerance = 1e-02))
> 
> ## eval_log_curve_exponential
> time <- seq.int(0, 100, by = 1)
> r <- log(2) / 20
> c0 <- 100
> ans <- getReport("eval_log_curve_exponential",
+                  time = time,
+                  log_r = log(r),
+                  log_c0 = log(c0))
> stopifnot(all.equal(ans, log(c0) + r * time))
> 
> ## eval_log_(curve|rt)_subexponential
> time <- seq.int(0, 100, by = 1)
> alpha <- log(2) / 20
> c0 <- 100
> p <- 0.95
> ans1 <- getReport("eval_log_curve_subexponential",
+                   time = time,
+                   log_alpha = log(alpha),
+                   log_c0 = log(c0),
+                   logit_p = qlogis(p))
> all.equal(ans1, log(c0) + log1p((1 - p) * alpha * time / c0^(1 - p)) / (1 - p))
[1] TRUE
> ans2 <- getReport("eval_log_rt_subexponential",
+                   log_curve = ans1,
+                   log_alpha = log(alpha),
+                   logit_p = qlogis(p))
> stopifnot(all.equal(ans2, log(alpha) - (1 - p) * ans1))
> 
> ## eval_log_(curve|rt)_gompertz
> time <- seq.int(0, 100, by = 1)
> alpha <- log(2) / 20
> tinfl <- 100
> K <- 25000
> ans1 <- getReport("eval_log_curve_gompertz",
+                   time = time,
+                   log_alpha = log(alpha),
+                   log_tinfl = log(tinfl),
+                   log_K = log(K))
> stopifnot(all.equal(ans1, log(K) - exp(-alpha * (time - tinfl))))
> ans2 <- getReport("eval_log_rt_gompertz",
+                   log_curve = ans1,
+                   log_alpha = log(alpha),
+                   log_K = log(K))
> stopifnot(all.equal(ans2, log(alpha) + log(log(K) - ans1)))
> 
> ## eval_log_(curve|rt)_logistic
> time <- seq.int(0, 100, by = 1)
> r <- log(2) / 20
> tinfl <- 100
> K <- 25000
> ans1 <- getReport("eval_log_curve_logistic",
+                   time = time,
+                   log_r = log(r),
+                   log_tinfl = log(tinfl),
+                   log_K = log(K))
> stopifnot(all.equal(ans1, log(K) - log1p(exp(-r * (time - tinfl)))))
> ans2 <- getReport("eval_log_rt_logistic",
+                   log_curve = ans1,
+                   log_r = log(r),
+                   log_K = log(K))
> stopifnot(all.equal(ans2, log(r) + log1p(-exp(ans1) / K)))
> 
> ## eval_log_(curve|rt)_richards
> time <- seq.int(0, 100, by = 1)
> r <- log(2) / 20
> tinfl <- 100
> K <- 25000
> a <- 1.005
> ans1 <- getReport("eval_log_curve_richards",
+                   time = time,
+                   log_r = log(r),
+                   log_tinfl = log(tinfl),
+                   log_K = log(K),
+                   log_a = log(a))
> stopifnot(all.equal(ans1, log(K) - log1p(a * exp(-a * r * (time - tinfl))) / a))
> ans2 <- getReport("eval_log_rt_richards",
+                   log_curve = ans1,
+                   log_r = log(r),
+                   log_K = log(K),
+                   log_a = log(a))
> stopifnot(all.equal(ans2, log(r) + log1p(-(exp(ans1) / K)^a)))
> 
> ## logspace_add_(baseline|offsets)
> time <- seq.int(0, 100, by = 1)
> log_curve <- log(100) + (log(2) / 20) * time
> log_diff_curve <- log(diff(exp(log_curve)))
> log_b <- log(2)
> log_w <- log(c(1, 0.9, 1.1, 0.8, 1.2, 0.7, 1.3))
> from <- 5L
> ans1 <- getReport("logspace_add_baseline",
+                   log_curve = log_curve,
+                   time = time,
+                   log_b = log_b)
> stopifnot(all.equal(ans1, log(exp(log_b) * time + exp(log_curve))))
> ans2 <- getReport("logspace_add_offsets",
+                   log_diff_curve = log_diff_curve,
+                   log_w = log_w,
+                   from = from)
> stopifnot(all.equal(ans2, log_diff_curve + rep_len(c(log_w[-seq_len(from)], log_w[seq_len(from)]), length(log_diff_curve))))
> 
> proc.time()
   user  system elapsed 
   5.67    0.23   96.81