context("mle")

test_that("mle works correctly for default Gaussian models", {

  AbaIdtmle <-mle(AbaloneIdt)
  
  BestM <- 1
  names(BestM) <- "NModCovC1"
  expect_equal(BestModel(AbaIdtmle),BestM)
  
  AllAbnames <- paste(names(AbaloneIdt),rep(c("MidP","LogR"),each=ncol(AbaloneIdt)),sep=".")

  Abmeans <- c( 0.5071875,    0.396145833,  0.166666667,  1.078052083,  0.45571875,   0.224302083,  0.3373125,
               -1.296114861, -1.534997069, -2.235509665,  0.095119706, -0.699579885, -1.43648703,  -1.176694026 )

  Abstddevs <- c( 0.126984892, 0.106306307, 0.099868838, 0.463170117, 0.204648509, 0.09817756, 0.159985436,
                  0.599929088, 0.646096586, 0.872978407, 1.098762777, 1.13577332, 1.108493279, 1.229327926 )
  
  Abcor <- matrix(
            c( 1, 0.993982764, 0.46722427, 0.943204961, 0.806954256, 0.877049526, 0.8988425, 0.197447713, 0.215699101, 0.254071722, 0.791917195, 0.767629178, 0.776267869, 0.717014941,
               0.993982764, 1, 0.460916165, 0.949745049, 0.801663444, 0.883022118, 0.914306533, 0.156517814, 0.176048065, 0.214254783, 0.75904035, 0.737417731, 0.74382086, 0.677761255,
               0.46722427, 0.460916165, 1, 0.507651246, 0.534390369, 0.492541963, 0.377370092, 0.206446863, 0.211721929, 0.621096109, 0.451146829, 0.457607674, 0.438228753, 0.374981242,
               0.943204961, 0.949745049, 0.507651246, 1, 0.911504486, 0.949135574, 0.923149776, 0.220628059, 0.237353382, 0.279750403, 0.76107623, 0.751149887, 0.746209619, 0.679318984,
               0.806954256, 0.801663444, 0.534390369, 0.911504486, 1, 0.914368882, 0.747872639, 0.444048709, 0.461436449, 0.494848721, 0.810320958, 0.832231663, 0.803650431, 0.723120244,
               0.877049526, 0.883022118, 0.492541963, 0.949135574, 0.914368882, 1, 0.825795256, 0.358369732, 0.37524005, 0.398271111, 0.794065249, 0.787335437, 0.807205016, 0.715729116,
               0.8988425, 0.914306533, 0.377370092, 0.923149776, 0.747872639, 0.825795256, 1, 0.051592822, 0.068467193, 0.096948743, 0.639028743, 0.609158995, 0.608926341, 0.627478887,
               0.197447713, 0.156517814, 0.206446863, 0.220628059, 0.444048709, 0.358369732, 0.051592822, 1, 0.992566651, 0.84448339, 0.696977234, 0.739845117, 0.717739093, 0.667198309,
               0.215699101, 0.176048065, 0.211721929, 0.237353382, 0.461436449, 0.37524005, 0.068467193, 0.992566651, 1, 0.835603923, 0.694992692, 0.741356203, 0.717499093, 0.665089517,
               0.254071722, 0.214254783, 0.621096109, 0.279750403, 0.494848721, 0.398271111, 0.096948743, 0.84448339, 0.835603923, 1, 0.692460681, 0.721881721, 0.702482351, 0.666151817,
               0.791917195, 0.75904035, 0.451146829, 0.76107623, 0.810320958, 0.794065249, 0.639028743, 0.696977234, 0.694992692, 0.692460681, 1, 0.990235553, 0.994727084, 0.960669307,
               0.767629178, 0.737417731, 0.457607674, 0.751149887, 0.832231663, 0.787335437, 0.609158995, 0.739845117, 0.741356203, 0.721881721, 0.990235553, 1, 0.986396398, 0.930830052,
               0.776267869, 0.74382086, 0.438228753, 0.746209619, 0.803650431, 0.807205016, 0.608926341, 0.717739093, 0.717499093, 0.702482351, 0.994727084, 0.986396398, 1, 0.952600258,
               0.717014941, 0.677761255, 0.374981242, 0.679318984, 0.723120244, 0.715729116, 0.627478887, 0.667198309, 0.665089517, 0.666151817, 0.960669307, 0.930830052, 0.952600258, 1),
            nrow=2*ncol(AbaloneIdt),ncol=2*ncol(AbaloneIdt)
  )

  names(Abmeans) <- names(Abstddevs) <- rownames(Abcor) <- colnames(Abcor) <- AllAbnames
  
  expect_equal(mean(AbaIdtmle),Abmeans)
  expect_equal(sd(AbaIdtmle),Abstddevs)
  expect_equal(cor(AbaIdtmle),Abcor)

} )


test_that("mle computes correct standar errors for default Gaussian models", {
  
  for (Cv in 1:4) {
    AbaIdtmle <-mle(AbaloneIdt, CovCase = Cv)

    n <-  nrow(AbaloneIdt) 
    q <-  2*ncol(AbaloneIdt)  # Total number of MidPoints and LogRanges 
  
    AbmeanStder <- sd(AbaIdtmle) / sqrt(n) 
    expect_equal(stdEr(AbaIdtmle)$mu,AbmeanStder)
    vcovb_AbmeanStder <- sqrt(diag(vcov(AbaIdtmle)[1:q,1:q]))
    names(vcovb_AbmeanStder) <- names(AbmeanStder) <- NULL
    expect_equal(vcovb_AbmeanStder,AbmeanStder)
  
    mlecov <- var(AbaIdtmle)
    mlecov[mlecov==0.] <- NA
    mlevar <- diag(mlecov)
    AbcovStder <- sqrt( (mlecov^2 + outer(mlevar,mlevar)) / (n-1) )   
    expect_equal(stdEr(AbaIdtmle)$Sigma,AbcovStder)

    if (Cv==1) {  # Implement later vcov checks for other configurations
      vcovb_AbcovStder <- matrix(nrow=q,ncol=q)
      vcovb_AbcovStder[lower.tri(vcovb_AbcovStder,diag=TRUE)] <- sqrt(diag(vcov(AbaIdtmle)[-(1:q),-(1:q)]))
      vcovb_AbcovStder[upper.tri(vcovb_AbcovStder)] <- t(vcovb_AbcovStder)[upper.tri(t(vcovb_AbcovStder))]
      dimnames(vcovb_AbcovStder) <- dimnames(AbcovStder)
      expect_equal(vcovb_AbcovStder,AbcovStder)
    }  
  }  
  
} )