R Under development (unstable) (2023-12-13 r85679 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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(aster)
> 
>  set.seed(42)
> 
>  pred <- c(0, 1, 1, 2)
>  fam <- c(1, 1, 3, 2)
>  famlist <- fam.default()
> 
>  nnode <- length(pred)
>  nind <- 7
> 
>  x <- rbind(c(1, 1, 2, 4),
+             c(1, 1, 3, 0),
+             c(1, 1, 5, 2),
+             c(1, 0, 1, 0),
+             c(1, 0, 4, 0),
+             c(0, 0, 0, 0),
+             c(0, 0, 0, 0))
> 
>  nrow(x) == nind
[1] TRUE
>  ncol(x) == nnode
[1] TRUE
> 
>  root <- matrix(1, nind, nnode)
> 
>  theta <- rnorm(nind * nnode, 0, .33)
>  theta <- matrix(theta, nind, nnode)
>  storage.mode(theta) <- "double"
> 
>  aster:::setfam(fam.default())
> 
>  out <- .C(aster:::C_aster_theta2phi,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      theta = theta,
+      phi = matrix(as.double(0), nind, nnode))
> 
>  my.phi <- theta
>  # poisson cumulant function
>  my.phi[ , 2] <- my.phi[ , 2] - exp(theta[ , 4])
>  # non-zero poisson cumulant function
>  my.phi[ , 1] <- my.phi[ , 1] - log(exp(exp(theta[ , 3])) - 1)
>  # bernoulli cumulant function
>  my.phi[ , 1] <- my.phi[ , 1] - log(1 + exp(theta[ , 2]))
> 
>  all.equal(out$phi, my.phi)
[1] TRUE
> 
>  tout <- .C(aster:::C_aster_phi2theta,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      phi = out$phi,
+      theta = matrix(as.double(0), nind, nnode))
> 
>  all.equal(tout$theta, theta)
[1] TRUE
> 
>  storage.mode(x) <- "double"
>  storage.mode(root) <- "double"
> 
>  cout <- .C(aster:::C_aster_theta2ctau,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      theta = theta,
+      ctau = matrix(as.double(0), nind, nnode))
> 
>  my.ctau <- theta
>  # poisson mean function
>  my.ctau[ , 4] <- exp(theta[ , 4])
>  # non-zero poisson mean function
>  foo <- exp(theta[ , 3])
>  my.ctau[ , 3] <- foo / (1 - exp(- foo))
>  # bernoulli cumulant function
>  my.ctau[ , 2] <- 1 / (1 + exp(- theta[ , 2]))
>  # bernoulli cumulant function
>  my.ctau[ , 1] <- 1 * 1 / (1 + exp(- theta[ , 1]))
> 
>  all.equal(cout$ctau, my.ctau)
[1] TRUE
> 
>  xout <- .C(aster:::C_aster_xpred,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      x = x,
+      root = root,
+      xpred = matrix(as.double(0), nind, nnode))
> 
>  my.xpred <- x
>  my.xpred[ , 1] <- root[ , 1]
>  my.xpred[ , 2] <- x[ , 1]
>  my.xpred[ , 3] <- x[ , 1]
>  my.xpred[ , 4] <- x[ , 2]
> 
>  all.equal(xout$xpred, my.xpred)
[1] TRUE
> 
>  tout <- .C(aster:::C_aster_ctau2tau,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      root = root,
+      ctau = cout$ctau,
+      tau = matrix(as.double(0), nind, nnode))
> 
>  my.tau <- my.ctau
>  my.tau[ , 1] <- my.xpred[ , 1] * my.ctau[ , 1]
>  my.tau[ , 2] <- my.tau[ , 1] * my.ctau[ , 2]
>  my.tau[ , 3] <- my.tau[ , 1] * my.ctau[ , 3]
>  my.tau[ , 4] <- my.tau[ , 2] * my.ctau[ , 4]
> 
>  all.equal(tout$tau, my.tau)
[1] TRUE
> 
>  ##########
> 
>  wout <- .C(aster:::C_aster_theta2whatsis,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      deriv = as.integer(0),
+      theta = as.double(theta),
+      result = matrix(as.double(0), nind, nnode))
> 
>  my.what <- NaN * wout$result
>  for (i in 1:nind)
+      for (j in 1:nnode) {
+          ifam <- fam[j]
+          my.what[i, j] <- famfun(famlist[[ifam]], 0, theta[i, j])
+      }
> 
>  all.equal(wout$result, my.what)
[1] TRUE
> 
>  ##########
> 
>  aster:::setfam(fam.default())
> 
>  wout <- .C(aster:::C_aster_theta2whatsis,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      deriv = as.integer(1),
+      theta = as.double(theta),
+      result = matrix(as.double(0), nind, nnode))
> 
>  my.what <- NaN * wout$result
>  for (i in 1:nind)
+      for (j in 1:nnode) {
+          ifam <- fam[j]
+          my.what[i, j] <- famfun(famlist[[ifam]], 1, theta[i, j])
+      }
>  all.equal(wout$result, my.what)
[1] TRUE
> 
>  ##########
> 
>  aster:::setfam(fam.default())
> 
>  wout <- .C(aster:::C_aster_theta2whatsis,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      deriv = as.integer(2),
+      theta = as.double(theta),
+      result = matrix(as.double(0), nind, nnode))
> 
>  my.what <- NaN * wout$result
>  for (i in 1:nind)
+      for (j in 1:nnode) {
+          ifam <- fam[j]
+          my.what[i, j] <- famfun(famlist[[ifam]], 2, theta[i, j])
+      }
> 
>  all.equal(wout$result, my.what)
[1] TRUE
> 
>  ##########
> 
>  aster:::setfam(fam.default())
> 
>  vout <- .C(aster:::C_aster_tt2var,
+      nind = as.integer(nind),
+      nnode = as.integer(nnode),
+      pred = as.integer(pred),
+      fam = as.integer(fam),
+      x = as.double(x),
+      root = as.double(root),
+      theta = as.double(theta),
+      tau = as.double(tout$tau),
+      var = matrix(as.double(0), nind, nnode))
> 
>  my.var <- NaN * vout$var
>  for (i in 1:nind)
+      for (j in 1:nnode) {
+          ifam <- fam[j]
+          thefam <- famlist[[ifam]]
+          k <- pred[j]
+          if (k > 0) {
+              my.var[i, j] <- famfun(thefam, 2, theta[i, j]) * tout$tau[i, k] +
+                  famfun(thefam, 1, theta[i, j])^2 * my.var[i, k]
+          } else {
+              my.var[i, j] <- famfun(thefam, 2, theta[i, j]) * xout$xpred[i, j]
+          }
+      }
>  all.equal(vout$var, my.var)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.17    0.06    0.21