generate.summaries <- function(result) {
  samples <- result$samples
  data <- as.matrix(samples)
  data <- data[, grep("^d\\.", colnames(data))]
  if (result$model$type == "consistency") {
    list(
      effectiveSize = effectiveSize(samples),
      summary       = summary(samples),
      cov           = cov(as.matrix(data)),
      ranks         = rank.probability(result),
      dic           = result$deviance[c('Dbar', 'pD', 'DIC', 'data points')]
    )
  } else {
    list(
      effectiveSize = effectiveSize(samples),
      summary       = summary(samples),
      cov           = cov(as.matrix(data)),
      dic           = result$deviance[c('Dbar', 'pD', 'DIC', 'data points')]
    )
  }
}

formatError <- function(name, s1, s2, test) {
  str.s1 <- paste(capture.output(print(s1)), collapse="\n")
  str.s2 <- paste(capture.output(print(s2)), collapse="\n")
  str.test <- paste(capture.output(print(test)), collapse="\n")

  paste(name,
        "=========== Expected:       ===========", str.s1,
        "=========== Actual:         ===========", str.s2,
        "=========== Test statistic: ===========", str.test, sep="\n\n")
}

compare.summaries <- function(s1, s2) {
  stopifnot(names(s1$effectiveSize) == names(s2$effectiveSize))
  d.idx <- grep("^d\\.", names(s1$effectiveSize))
  sd.idx <- grep("^sd\\.", names(s1$effectiveSize))

  # Force statistics to be matrix
  if (is.vector(s1$summary$statistics)) s1$summary$statistics <- t(s1$summary$statistics)
  if (is.vector(s2$summary$statistics)) s2$summary$statistics <- t(s2$summary$statistics)

  # Test equality of means
  mu1 <- s1$summary$statistics[, 'Mean']
  mu2 <- s2$summary$statistics[, 'Mean']
  se1 <- s1$summary$statistics[, 'Time-series SE']
  se2 <- s2$summary$statistics[, 'Time-series SE']
  test <- pnorm(mu1 - mu2, 0, sqrt(se1^2 + se2^2))
  expect(all(test > 0.025), formatError("Means were not equal", mu1, mu2, test))

  # Test equality of variance
  sd1 <- s1$summary$statistics[d.idx, 'SD']
  sd2 <- s2$summary$statistics[d.idx, 'SD']
  en1 <- s1$effectiveSize[d.idx]
  en2 <- s2$effectiveSize[d.idx]
  test <- pf(sd1^2 / sd2^2, en1, en2)
  expect(all(test > 0.025), formatError("Variances were not equal", sd1^2, sd2^2, test))

  # TODO: multivariate test for equality of means
  # TODO: compare covariance matrices

  # TODO: compare quantiles for the standard deviation using http://www.jstor.org/stable/2673594

  # Test equality of rank probabilities
  if (!is.null(s1$ranks)) {
    thin <- s2$summary$thin
    n.adapt <- s2$summary$start - thin
    n.iter <- s2$summary$end - n.adapt
    test <- sapply(rownames(s2$ranks), function(alt) {
      n.sample <- min(s2$effectiveSize[d.idx])
      x <- round(s2$ranks[alt, ] * n.sample)
      p <- s1$ranks[alt, ]

      # Continuity correction because too many zero cells cause the test to fail
      x <- x + 1
      p <- (p + 1/n.sample) / (1 + nrow(s2$ranks)/n.sample)

      test <- chisq.test(x, p=p, rescale.p=TRUE, simulate.p.value=TRUE)
      c('statistic'=unname(test$statistic), 'p.value'=test$p.value)
    })
    expect(all(test['p.value',] > 0.025), formatError("Rank probabilities were not equal", s1$ranks, s2$ranks, test))
  }

  # Test equality of deviance statistics
  if (!is.null(s1$dic)) {
    n <- min(s1$effectiveSize, s2$effectiveSize)
    expect_equal(s1$dic[['data points']], s2$dic[['data points']])
    # deviance should follow an approximate Chi-squared distribution with variance 2*df
    se <- sqrt(2 * s1$dic[['data points']] / n)
    v1 <- unlist(s1$dic[c('Dbar', 'pD')])
    v2 <- unlist(s2$dic[c('Dbar', 'pD')])
    test <- pnorm(v1 - v2, 0, se)
    expect(all(test > 0.025), formatError("Model fit statistics were not equal", v1, v2, test))
  }
}

replicate.example <- function(name, network, type="consistency", linearModel="random", likelihood="binom", link="logit") {
  s1 <- dget(paste0("../data/", name, '.summaries.txt'))

  n.chain <- s1$summary$nchain
  thin <- s1$summary$thin
  n.adapt <- s1$summary$start - thin
  n.iter <- s1$summary$end - n.adapt

  model <- if (exists("powerAdjustMode")) {
    network$studies <- data.frame(study=mtc.studies.list(network)[['values']], x=1)
    mtc.model(network, type=type,
      likelihood=likelihood, link=link,
      linearModel=linearModel,
      n.chain=4,
      powerAdjust='x')
  } else {
    mtc.model(network, type=type,
      likelihood=likelihood, link=link,
      linearModel=linearModel,
      n.chain=4)
  }
  capture.output(
    result <- mtc.run(model, n.adapt=n.adapt, n.iter=n.iter, thin=thin)
  )
  s2 <- generate.summaries(result)
  list(s1=s1, s2=s2)
}

verify.example <- function(example, sampler) {
  x <- replicate.example(example, sampler)
  compare.summaries(x$s1, x$s2)
}