#devtools::test("dae")
context("Ameasures")

cat("#### Test for mat.Vpred and designAmeasures using Smith's small example\n")
test_that("SmithSmall", {
  skip_on_cran()
  library(dae)
  #'# Script to investigate the reduced design from Smith et al. (2015)

  #'## Set up designs
  mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
  field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
  field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                              levels = c("Y","W","G","M","D","E"))
  start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
  rownames(start.design) <- NULL
  start.design

  #'## Set gammas
  terms <- c("Variety", "Frep", "Frep:Fplot", "Mrep", "Mrep:Mday", "Mrep:Mday:Mord")
  gammas <- c(1, 0.1, 0.2, 0.3, 0.2, 1)
  names(gammas) <- terms
  print(gammas)
  
  #'## Examine A-optimality measures
  #'### Set up functions
  "invInf" <- function(design, gammas, ...)
  {
    #'## Set up matrices
    W <- model.matrix(~ -1 + Variety, design)
    ZMr <- model.matrix(~ -1 + Mrep, design)
    ZMrMd <- model.matrix(~ -1 + Mrep:Mday, design)
    ZFr <- model.matrix(~ -1 + Frep, design)
    ZFrFp <- model.matrix(~ -1 + Frep:Fplot, design)
    t <- length(levels(design$Variety))
    ng <- ncol(W)
    Vu <- gammas["Mrep"] * ZMr %*% t(ZMr) + gammas["Mrep:Mday"] * ZMrMd %*% t(ZMrMd) + 
      gammas["Frep"] * ZFr %*% t(ZFr) + gammas["Frep:Fplot"] * ZFrFp %*% t(ZFrFp)
    R <- diag(1, nrow(design))
    
    #'# Calculate variance matrix
    Vp <- mat.Vpred(W = W, Vu=Vu, R=R, ...)
    
    return(Vp)
  }
  
  
  #'## Calculate A measures
  Vp <- invInf(start.design, gammas)
  testthat::expect_true(all(abs(designAmeasures(Vp) - 1.521905) < 1e-6))
  n <- nrow(start.design)
  Vp.reduc <- invInf(start.design, gammas, 
                     eliminate = projector(matrix(1, nrow = n, ncol = n)/n))
  testthat::expect_true(all(abs(rowSums(Vp.reduc)) < 1e-10))
  testthat::expect_true(abs(designAmeasures(Vp.reduc) - 
                              sum(diag(Vp.reduc)) * 2 / (nrow(Vp.reduc) - 1)) < 1.e-10)
  
  #'## Investigate the anatomy
  start2.canon <- designAnatomy(list(lab = ~ Mrep/Mday/Mord, 
                                     plot = ~ Frep/Fplot),
                                data = start.design, omit.projectors = "pcanon")
  testthat::expect_equal(degfree(start2.canon$Q[[2]]$`Mord[Mrep:Mday]&Fplot[Frep]`), 5)
})  

cat("#### Test for Ameasures using Cochran&Cox PBIBD2\n")
test_that("PBIBD2Ameasures", {
  skip_on_cran()
  library(dae)
  #'# PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 2nd edn Wiley, New York"
  
  #'## Input the design and randomize"
  Treatments <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
  PBIBD2.lay <- designRandomize(allocated = Treatments,
                                recipient = list(Blocks =6, Units = 4),
                                nested.recipients = list(Units = "Blocks"),
                                seed = 98177)
  
  #'## Compute the anatomy
  PBIBD2.canon <- designAnatomy(formulae = list(unit = ~Blocks/Units,
                                                trt = ~ Treatments),
                                which.criteria = c('aeff', 'xeff', 'eeff','order'),
                                data = PBIBD2.lay, omit.projectors = "combined")
  
  
  testthat::expect_equal(PBIBD2.canon$Q[[1]]$Blocks$Treatments$adjusted$aefficiency, 0.25)
  testthat::expect_lt(abs(PBIBD2.canon$Q[[1]]$`Units[Blocks]`$Treatments$adjusted$aefficiency - 0.8824), 1e-04)
  
  #'## Get Ameasure
  #'## Set up matrices
  W <- model.matrix(~ -1 + Treatments, PBIBD2.lay)
  Vu <- with(PBIBD2.lay, fac.vcmat(Blocks, 5))
  R <- diag(1, nrow(PBIBD2.lay))
  
  #'## Calculate the variance matrix of the predictions
  Vp <- mat.Vpred(W = W, Vu=Vu, R=R)
  testthat::expect_true(abs(designAmeasures(Vp) - 0.5625) < 1e-6)

  #Test when no random terms
  Vpr <- mat.Vpredicts(target = ~ -1 + Treatments, fixed = ~ Blocks -1, design = PBIBD2.lay)
  testthat::expect_true(abs(designAmeasures(Vpr) - 0.5666667) < 1e-6)
  
    
  #'## Get reduced variance matix of predictions
  unit.str <- pstructure(~Blocks/Units, grandMean = TRUE, data = PBIBD2.lay)
  Q.B <- projector(unit.str$Q$Mean + unit.str$Q$Blocks)
  testthat::expect_equal(degfree(Q.B), 6)
  Vp.reduc <- mat.Vpred(W = W, Vu=Vu, R=R, eliminate = Q.B)
  testthat::expect_true(all(abs(rowSums(Vp.reduc)) < 1e-10))
  testthat::expect_true(abs(designAmeasures(Vp.reduc) - 
                              2/(4*PBIBD2.canon$Q[[1]]$`Units[Blocks]`$Treatments$adjusted$aefficiency)) 
                        < 1e-10)
})

cat("#### Test for mat.Vpredicts using using Smith's small example\n")
test_that("SmithSmallVpredicts", {
  skip_on_cran()
  library(dae)
  #'# Script to investigate the reduced design from Smith et al. (2015)
  
  #'## Set up designs
  mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
  field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
  field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                              levels = c("Y","W","G","M","D","E"))
  start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
  rownames(start.design) <- NULL
  start.design
  
  #'## Set gammas
  terms <- c("Variety", "Frep", "Frep:Fplot", "Mrep", "Mrep:Mday", "Mrep:Mday:Mord")
  gammas <- c(1, 0.1, 0.2, 0.3, 0.2, 1)
  names(gammas) <- terms
  print(gammas)
  
  #'## Set up matrices for fixed Variety effects
  W <- model.matrix(~ -1 + Variety, start.design)
  ZMr <- model.matrix(~ -1 + Mrep, start.design)
  ZMrMd <- model.matrix(~ -1 + Mrep:Mday, start.design)
  ZFr <- model.matrix(~ -1 + Frep, start.design)
  ZFrFp <- model.matrix(~ -1 + Frep:Fplot, start.design)
  t <- length(levels(start.design$Variety))
  ng <- ncol(W)
  Vu <- gammas["Mrep"] * ZMr %*% t(ZMr) + gammas["Mrep:Mday"] * ZMrMd %*% t(ZMrMd) + 
    gammas["Frep"] * ZFr %*% t(ZFr) + gammas["Frep:Fplot"] * ZFrFp %*% t(ZFrFp)
  R <- diag(1, nrow(start.design))
  
  #'# Calculate variance matrix
  Vp <- mat.Vpred(W = W, Vu=Vu, R=R)
 
  ## Specify matrices to calculate the variance matrix of the predicted fixed Variety effects 
  W <- model.matrix(~ -1 + Variety, start.design)
  Vu <- with(start.design, fac.vcmat(Mrep, gammas["Mrep"]) + 
               fac.vcmat(fac.combine(list(Mrep,Mday)), gammas["Mrep:Mday"]) + 
               fac.vcmat(Frep, gammas["Frep"]) + 
               fac.vcmat(fac.combine(list(Frep,Fplot)), gammas["Frep:Fplot"]))
  R <- diag(1, nrow(start.design))
  
  ## Calculate variance matrix
  Vp <- mat.Vpredicts(target = W, random = Vu, R = R, design = start.design)
  testthat::expect_equal(attr(Vp, which = "rank"), 5)
             
  ## Calculate information matrix
  C <- mat.Vpredicts(target = W, random = Vu, R = R, result = "inform", design = start.design)
  testthat::expect_equal(attr(C, which = "rank"), 5)
  
  #'## Calculate A measures
  AVpred <- designAmeasures(Vp)
  testthat::expect_true(all(abs(AVpred - 1.521905) < 1e-6))
  
  Vp <- mat.Vpredicts(target = W, random = Vu, R = R, design = start.design)
  AVpredict <- designAmeasures(Vp)
  testthat::expect_true(all(abs(AVpredict - AVpred) < 1e-6))
  
  Vp <- mat.Vpredicts(target = ~ -1 + Variety, 
                      fixed = ~ 1, 
                      random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                      G = as.list(gammas[c(4,5,2,3)]), 
                      R = R, design = start.design)
  Aform <- designAmeasures(Vp)
  testthat::expect_true(all(abs(Aform - AVpred) < 1e-6))
  
  ## Calculate the variance matrix of the predicted random Variety effects using formulae
  Vp <- mat.Vpredicts(target = ~ -1 + Variety, Gt = 1, 
                      fixed = ~ 1, 
                      random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                      G = as.list(gammas[c(4,5,2,3)]), 
                      R = R, design = start.design)
  Aran <- designAmeasures(Vp)
  testthat::expect_true(all(abs(Aran - 0.8564816) < 1e-6))
  
  
  ## Calculate the variance matrix of the predicted fixed Variety effects, 
  ## elminating the grand mean
  n <- nrow(start.design)
  Vp.reduc <- mat.Vpredicts(target = ~ -1 + Variety, 
                            random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                            G = as.list(gammas[c(4,5,2,3)]), 
                            eliminate = projector(matrix(1, nrow = n, ncol = n)/n), 
                            design = start.design)
  testthat::expect_true(all(abs(rowSums(Vp.reduc)) < 1e-10))
  testthat::expect_true(abs(designAmeasures(Vp.reduc) - 
                              sum(diag(Vp.reduc)) * 2 / (nrow(Vp.reduc) - 1)) < 1.e-10)
})  


cat("#### Test for Ameasures using Cochran&Cox PBIBD2\n")
test_that("PBIBD2Ameasures", {
  skip_on_cran()
  library(dae)
  
  #Input systematic design
  RIBDmultir.sys <- cbind(fac.gen(list(Blocks = 8, Units = 6)),
                          Treats = factor(c(3,4,5,6,1,2, 4,2,6,3,5,1,
                                            5,3,2,1,4,2, 6,1,3,2,2,4,
                                            3,5,1,2,1,4, 1,2,2,4,3,3,
                                            2,4,1,5,3,6, 2,3,4,1,2,5)))
  
  #Randomize the design
  RIBDmultir.lay <- designRandomize(recipient = RIBDmultir.sys[c("Blocks", "Units")],
                                    nested.recipients = list(Units = "Blocks"), 
                                    allocated = RIBDmultir.sys["Treats"], 
                                    seed = 71571)

  #Check the properties of the desing
  RIBDmultir.canon <- designAnatomy(list(unit = ~ Blocks/Units, 
                                         trts = ~ Treats), 
                                    data = RIBDmultir.lay)
  summary(RIBDmultir.canon)

  #Get the canonical efficiency factors and the replication vector
  effs <- efficiencies(RIBDmultir.canon)
  effs <- effs[[1]]$`Units[Blocks]`$Treats
  testthat::expect_true(all(abs(effs - c(1.0000000,0.9907854,0.9814815,0.9408367,0.9017113)) < 1e-05))
  r <- as.vector(table(RIBDmultir.lay$Treats))
  rmean <- mean(r)
  testthat::expect_equal(rmean, 8)
  
  #Calculate A-measures using the canonical efficiency factors
  e_a <- harmonic.mean(effs)
  testthat::expect_true(abs(e_a - 0.9615284) < 1e-05)
  
  #Get predictions variance matrix
  Vpred <- mat.Vpredicts(target = ~ -1 + Treats,
                         fixed = ~ 1 + Blocks, 
                         design = RIBDmultir.lay)
  
  #Compute A-measures from predictions variance matrix
  apv <- designAmeasures(Vpred)
  g_a <- 2/apv/rmean
  f_a <- g_a * (rmean/harmonic.mean((r)))
  testthat::expect_true(abs(f_a - 0.9528458) < 1e-05)
    e_a <- 2/designAmeasures(Vpred, replications = r)/rmean
  testthat::expect_true(abs(e_a - 0.9615284) < 1e-05)
  
})