test_that("mdata", {
  local_bru_safe_inla()

  sim.poisson <- function(prob, m) {
    stopifnot(length(prob) == length(m) && length(prob) > 0)
    n <- length(m)
    y <- numeric(n)
    event <- (runif(n) < prob)
    idx.zero <- which(event)
    idx.non.zero <- which(!event)
    y[idx.zero] <- 0
    y[idx.non.zero] <- rpois(length(idx.non.zero), lambda = m[idx.non.zero])
    return(y)
  }
  ## chose link-function to use for the zero-inflation probability
  link.simple <- "logit"
  inv.link <- INLA::inla.link.invlogit
  ## link.simple <- "probit"
  ## inv.link <- inla.link.invprobit
  ## link.simple <- "cloglog"
  ## inv.link <- inla.link.invcloglog
  n <- 1000
  z <- rnorm(n, sd = 0.3)
  x <- rnorm(n, sd = 0.2)
  xx <- rnorm(n, sd = 0.3)
  zz <- rnorm(n, sd = 0.2)
  E <- runif(n, min = 0.8, max = 1 / 0.8)
  beta <- c(1, 1.1, 2.1, 0, -2, 1.2, 2.2, 0)
  eta2 <- beta[1] + beta[2] * xx + beta[3] * zz + beta[4] * xx * zz
  eta1 <- beta[5] + beta[6] * x + beta[7] * z + beta[8] * x * z
  prob <- inv.link(eta1)
  m <- E * exp(eta2)
  ok <- FALSE
  while (!ok) {
    y <- sim.poisson(prob, m)
    ok <- !all(y == 0)
  }
  ## head(data.frame(y, E, x, z, xx, zz))
  suppressWarnings(
    r <- INLA::inla(
      INLA::inla.mdata(cbind(y, E), cbind(1, x, z, x * z)) ~ 1 + xx + zz + xx * zz,
      family = "0poisson",
      data = data.frame(y, E, x, z, xx, zz),
      control.fixed = list(prec = 1, prec.intercept = 1),
      control.compute = list(cpo = TRUE),
      control.family = list(
        link.simple = link.simple,
        hyper = list(
          beta1 = list(param = c(0, 1)),
          beta2 = list(param = c(0, 1)),
          beta3 = list(param = c(0, 1)),
          beta4 = list(param = c(0, 1)),
          beta5 = list(param = c(0, 1))
        )
      )
    )
  )
  #  y2 <- INLA::inla.mdata(cbind(y, -1, E/2),
  #                         cbind(1, x, z, x*z))
  #  r2 <- INLA::inla(
  #    y2 ~ 1 + xx + zz + xx*zz,
  #    family = "0poisson",
  #    data = data.frame(y, E, x, z, xx, zz),
  #    control.fixed = list(prec = 1, prec.intercept = 1),
  #    control.compute = list(cpo = TRUE),
  #    control.family = list(
  #      list(link.simple = link.simple,
  #                          hyper = list(beta1 = list(param = c(0, 1)),
  #                                       beta2 = list(param = c(0, 1)),
  #                                       beta3 = list(param = c(0, 1)),
  #                                       beta4 = list(param = c(0, 1)),
  #                                       beta5 = list(param = c(0, 1)))))
  #  )
  #

  if (FALSE) {
    r_bru <- bru(
      ~ 0 + fix(~ 1 + xx + zz + xx * zz, model = "fixed", prec.linear = 1),
      like(
        INLA::inla.mdata(cbind(y, E), cbind(1, x, z, x * z)) ~ .,
        family = "0poisson",
        data = data.frame(y, E, x, z, xx, zz),
        control.family = list(
          link.simple = link.simple,
          hyper = list(
            beta1 = list(param = c(0, 1)),
            beta2 = list(param = c(0, 1)),
            beta3 = list(param = c(0, 1)),
            beta4 = list(param = c(0, 1)),
            beta5 = list(param = c(0, 1))
          )
        )
      ),
      options = list(
        control.compute = list(cpo = TRUE)
      )
    )
  }


  if (FALSE) {
    rr <- inla(
      inla.mdata(cbind(y, E), cbind(1, xx, zz, xx * zz)) ~ 1 + x + z + x * z,
      family = "0poissonS",
      data = data.frame(y, E, x, z, xx, zz),
      control.fixed = list(prec = 1, prec.intercept = 1),
      control.compute = list(cpo = TRUE),
      ## in this case we need to define link.simple as the main link
      control.family = list(
        control.link = list(model = link.simple),
        hyper = list(
          beta1 = list(param = c(0, 1)),
          beta2 = list(param = c(0, 1)),
          beta3 = list(param = c(0, 1)),
          beta4 = list(param = c(0, 1)),
          beta5 = list(param = c(0, 1))
        )
      )
    )
    summary(r)
    summary(rr)
    res <- cbind(
      "beta" = beta,
      "0poisson" = c(r$summary.fixed$mean, r$summary.hyperpar$mean),
      "0poissonS" = c(rr$summary.hyperpar$mean, rr$summary.fixed$mean)
    )
    res <- cbind(res,
      diff = (res[, 2] - beta),
      diffS = (res[, 3] - beta),
      "diff/sd" = (res[, 2] - beta) / c(r$summary.fixed$sd, r$summary.hyperpar$sd),
      "diffS/sd" = (res[, 3] - beta) / c(rr$summary.hyperpar$sd, rr$summary.fixed$sd)
    )
    mm <- nrow(res) %/% 2
    rownames(res) <- c(paste0("beta", 1:mm, ".poisson"), paste0("beta", 1:mm, ".prob"))
    print(round(dig = 2, res))
  }
})