tests/testthat/test-cholperm.R

In nlmixr2/rxode2: Facilities for Simulating from ODE-Based Models

rxTest({
  # These functions are taken from TruncatedNormal for testing
  .rx <- loadNamespace("rxode2")

  lnNpr <-
    function(a, b) { ## computes ln(P(a<Z<b))
      ## where Z~N(0,1) very accurately for any 'a', 'b'
      p <- rep(0, length(a))
      ## case b>a>0
      I <- a > 0
      if (any(I)) {
        pa <- pnorm(a[I], lower.tail = FALSE, log.p = TRUE)
        pb <- pnorm(b[I], lower.tail = FALSE, log.p = TRUE)
        p[I] <- pa + log1p(-exp(pb - pa))
      }
      ## case a<b<0
      idx <- b < 0
      if (any(idx)) {
        pa <- pnorm(a[idx], log.p = TRUE)
        pb <- pnorm(b[idx], log.p = TRUE)
        p[idx] <- pb + log1p(-exp(pa - pb))
      }
      ## case a<0<b
      I <- !I & !idx
      if (any(I)) {
        pa <- pnorm(a[I])
        pb <- pnorm(b[I], lower.tail = FALSE)
        p[I] <- log1p(-pa - pb)
      }
      return(p)
    }

  cholperm <-
    function(Sig, l, u) {
      ##  Computes permuted lower Cholesky factor L for Sig
      ##  by permuting integration limit vectors l and u.
      ## Outputs perm, such that Sig(perm,perm)=L%*%t(L).
      ##
      ## Reference:
      ##  Gibson G. J., Glasbey C. A., Elston D. A. (1994),
      ##  "Monte Carlo evaluation of multivariate normal integrals and
      ##  sensitivity to variate ordering",
      ##  In: Advances in Numerical Methods and Applications, pages 120--126
      eps <- 10^-10 # round-off error tolerance
      d <- length(l)
      perm <- 1:d # keep track of permutation
      L <- matrix(0, d, d)
      z <- rep(0, d)
      for (j in 1:d) {
        pr <- rep(Inf, d) # compute marginal prob.
        I <- j:d # search remaining dimensions
        D <- diag(Sig)
        if (j > 2) {
          s <- D[I] - L[I, 1:(j - 1)]^2 %*% rep(1, j - 1)
        } else if (j == 2) {
          s <- D[I] - L[I, 1]^2
        } else {
          s <- D[I]
        }
        s[s < 0] <- eps
        s <- sqrt(s)
        if (j > 2) {
          cols <- L[I, 1:(j - 1)] %*% z[1:(j - 1)]
        } else if (j == 2) {
          cols <- L[I, 1] * z[1]
        } else {
          cols <- 0
        }
        tl <- (l[I] - cols) / s
        tu <- (u[I] - cols) / s
        pr[I] <- lnNpr(tl, tu)
        # find smallest marginal dimension
        k <- which.min(pr)
        # flip dimensions k-->j
        jk <- c(j, k)
        kj <- c(k, j)
        Sig[jk, ] <- Sig[kj, ]
        Sig[, jk] <- Sig[, kj] # update rows and cols of Sig
        L[jk, ] <- L[kj, ] # update only rows of L
        l[jk] <- l[kj]
        u[jk] <- u[kj] # update integration limits
        perm[jk] <- perm[kj] # keep track of permutation
        # construct L sequentially via Cholesky computation
        s <- Sig[j, j] - sum(L[j, 1:(j - 1)]^2)
        if (s < (-0.001)) {
          stop("Sigma is not positive semi-definite")
        }
        s[s < 0] <- eps
        L[j, j] <- sqrt(s)
        if (j < d) {
          if (j > 2) {
            L[(j + 1):d, j] <- (Sig[(j + 1):d, j] - L[(j + 1):d, 1:(j - 1)] %*% L[j, 1:(j - 1)]) / L[j, j]
          } else if (j == 2) {
            L[(j + 1):d, j] <- (Sig[(j + 1):d, j] - L[(j + 1):d, 1] * L[j, 1]) / L[j, j]
          } else if (j == 1) {
            L[(j + 1):d, j] <- Sig[(j + 1):d, j] / L[j, j]
          }
        }
        ## find mean value, z(j), of truncated normal:
        tl <- (l[j] - L[j, 1:j] %*% z[1:j]) / L[j, j]
        tu <- (u[j] - L[j, 1:j] %*% z[1:j]) / L[j, j]
        w <- lnNpr(tl, tu) # aids in computing expected value of trunc. normal
        z[j] <- (exp(-.5 * tl^2 - w) - exp(-.5 * tu^2 - w)) / sqrt(2 * pi)
      }
      return(list(L = L, l = l, u = u, perm = perm))
    }

  gradpsi <-
    function(y, L, l, u) { # implements grad_psi(x) to find optimal exponential twisting;
      # assume scaled 'L' with zero diagonal;
      d <- length(u)
      c <- rep(0, d)
      x <- c
      mu <- c
      x[1:(d - 1)] <- y[1:(d - 1)]
      mu[1:(d - 1)] <- y[d:(2 * d - 2)]
      # compute now ~l and ~u
      c[-1] <- L[-1, ] %*% x
      lt <- l - mu - c
      ut <- u - mu - c
      # compute gradients avoiding catastrophic cancellation
      w <- lnNpr(lt, ut)
      pl <- exp(-0.5 * lt^2 - w) / sqrt(2 * pi)
      pu <- exp(-0.5 * ut^2 - w) / sqrt(2 * pi)
      P <- pl - pu
      # output the gradient
      dfdx <- -mu[-d] + t(t(P) %*% L[, -d])
      dfdm <- mu - x + P
      grad <- c(dfdx, dfdm[-d])
      # here compute Jacobian matrix
      lt[is.infinite(lt)] <- 0
      ut[is.infinite(ut)] <- 0
      dP <- (-P^2) + lt * pl - ut * pu # dPdm
      DL <- rep(dP, 1, d) * L
      mx <- -diag(d) + DL
      xx <- t(L) %*% DL
      mx <- mx[1:(d - 1), 1:(d - 1)]
      xx <- xx[1:(d - 1), 1:(d - 1)]
      if (d > 2) {
        Jac <- rbind(cbind(xx, t(mx)), cbind(mx, diag(1 + dP[1:(d - 1)])))
      } else {
        Jac <- rbind(cbind(xx, t(mx)), cbind(mx, 1 + dP[1:(d - 1)]))
        dimnames(Jac) <- NULL
      }
      f <- list(grad = grad, Jac = Jac)
    }

  nleq <-
    function(l, u, L) {
      d <- length(l)
      x <- rep(0, 2 * d - 2) # initial point for Newton iteration
      err <- Inf
      iter <- 0
      while (err > 10^-10) {
        f <- gradpsi(x, L, l, u)
        Jac <- f$Jac
        grad <- f$grad
        del <- solve(Jac, -grad) # Newton correction
        x <- x + del
        err <- sum(grad^2)
        iter <- iter + 1
        if (iter > 100) {
          stop("Covariance matrix is ill-conditioned and method failed")
        }
      }
      return(x)
    }

  test_that("cholperm", {

    rxWithSeed(12, {
      d <- 5

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)

      r1 <- cholperm(mcov, -(1:5), 1:5)
      r2 <- .rx$rxCholperm(mcov, -(1:5), 1:5)

      expect_equal(r1$L, r2$L)
      expect_equal(r1$l, r2$l)
      expect_equal(r1$u, r2$u)
      expect_equal(r1$perm, r2$perm + 1)

      r1 <- cholperm(mcov, 1:5, 2 * (1:5))
      r2 <- .rx$rxCholperm(mcov, 1:5, 2 * 1:5)

      expect_equal(r1$L, r2$L)
      expect_equal(r1$l, r2$l)
      expect_equal(r1$u, r2$u)
      expect_equal(r1$perm, r2$perm + 1)

      r1 <- cholperm(mcov, -2 * (1:5), -(1:5))
      r2 <- .rx$rxCholperm(mcov, -2 * (1:5), -(1:5))

      expect_equal(r1$L, r2$L)
      expect_equal(r1$l, r2$l)
      expect_equal(r1$u, r2$u)
      expect_equal(r1$perm, r2$perm + 1)

      ## microbenchmark::microbenchmark(cholperm(mcov, -2 * (1:5), -(1:5)), rxCholperm(mcov, -2 * (1:5), -(1:5)))
      ## microbenchmark::microbenchmark(microbenchmark::cholperm(mcov, -2 * (1:5), -(1:5)), rxCholperm(mcov, -2 * (1:5), -(1:5)))
    })

  })

  test_that("gradpsi", {

    rxWithSeed(12, {
      d <- 5

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)
      r2 <- .rx$rxCholperm(mcov, -(1:5), 1:5)

      d <- length(r2$l)
      y <- seq(0, 2 * d - 2)
      l <- r2$l
      u <- r2$u
      L <- r2$L

      r1 <- gradpsi(y, L, l, u)
      r2 <- .rx$rxGradpsi(y, L, l, u)

      expect_equal(r1$Jac, r2$Jac)
      expect_equal(r1$grad, r2$grad)

      r1 <- gradpsi(rep(1, 2 * d - 2), L, l, u)
      r2 <- .rx$rxGradpsi(rep(1, 2 * d - 2), L, l, u)

      expect_equal(r1$Jac, r2$Jac)
      expect_equal(r1$grad, r2$grad)

      r1 <- gradpsi(rep(-1, 2 * d - 2), L, l, u)
      r2 <- .rx$rxGradpsi(rep(-1, 2 * d - 2), L, l, u)

      expect_equal(r1$Jac, r2$Jac)
      expect_equal(r1$grad, r2$grad)

      ## microbenchmark::microbenchmark(gradpsi(rep(-1,2*d-2), L, l, u), .rx$rxGradpsi(rep(-1,2*d-2), L, l, u));

      d <- 2

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)
      r2 <- .rx$rxCholperm(mcov, -(1:d), 1:d)

      d <- length(r2$l)
      y <- seq(0, 2 * d - 2)
      l <- r2$l
      u <- r2$u
      L <- r2$L


      r1 <- gradpsi(rep(-1, 2 * d - 2), L, l, u)
      r2 <- .rx$rxGradpsi(rep(-1, 2 * d - 2), L, l, u)

      expect_equal(r1$Jac, r2$Jac)
      expect_equal(r1$grad, r2$grad)
    })

  })


  test_that("nleq", {

    rxWithSeed(12, {
      d <- 5

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)
      r2 <- .rx$rxCholperm(mcov, -3 * (1:d), 2 * (1:d))

      expect_equal(.rx$rxNleq(r2$l, r2$u, r2$L), nleq(r2$l, r2$u, r2$L))
      ## microbenchmark::microbenchmark(rxNleq(r2$l, r2$u, r2$L), nleq(r2$l, r2$u, r2$L))

      d <- 2

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)
      r2 <- .rx$rxCholperm(mcov, -2 * (1:d), 3 * 1:d)

      expect_equal(.rx$rxNleq(r2$l, r2$u, r2$L), nleq(r2$l, r2$u, r2$L))
    })

  })

  test_that("rxMvnrnd", {

    rxWithSeed(12, {
      d <- 5

      mu <- 1:d

      ## Creating covariance matrix
      tmp <- matrix(rnorm(d^2), d, d)
      mcov <- tcrossprod(tmp, tmp)

      out <- .rx$rxCholperm(mcov, -2 * (1:5), 1:5)

      Lfull <- out$L
      l <- out$l
      u <- out$u
      D <- diag(Lfull)
      perm <- out$perm
      if (any(D < 10^-10)) {
        warning("Method may fail as covariance matrix is singular!")
      }
      L <- Lfull / D
      u <- u / D
      l <- l / D # rescale
      L <- L - diag(d) # remove diagonal
      # find optimal tilting parameter via non-linear equation solver
      xmu <- nleq(l, u, L) # nonlinear equation solver
      x <- xmu[1:(d - 1)]
      mu <- xmu[d:(2 * d - 2)] # assign saddlepoint x* and mu*

      fun <- function(n) {
        r1 <- .rx$rxMvnrnd(5, L, l, u, mu)
        expect_equal(length(r1$logpr), 5)
        expect_true(all(!duplicated(r1$logpr)))
        expect_equal(length(r1$Z[1, ]), 5)
        expect_true(all(!duplicated(r1$Z[1, ])))
      }
      fun(2)
      fun(5)
      fun(10)
    })

  })
})