library(lme4)
library(testthat)
library(lattice)
testLevel <- if (nzchar(s <- Sys.getenv("LME4_TEST_LEVEL"))) as.numeric(s) else 1

options(nwarnings = 5000)# instead of 50, and then use  summary(warnings())

if (testLevel>1) {
### __ was ./profile_plots.R ___
fm1 <- lmer(Reaction~ Days + (Days|Subject), sleepstudy)
pfile <- system.file("testdata","tprfm1.RData", package="lme4")
if(file.exists(pfile)) print(load(pfile)) else withAutoprint({
    system.time( tpr.fm1 <- profile(fm1, optimizer="Nelder_Mead") ) ## 5 sec (2018); >= 50 warnings !?
    save(tpr.fm1, file= "../../inst/testdata/tprfm1.RData")
})
oo <- options(warn = 2) # {warnings are errors from here on}

if(!dev.interactive(orNone=TRUE)) pdf("profile_plots.pdf")
xyplot(tpr.fm1)
splom(tpr.fm1)
densityplot(tpr.fm1, main="densityplot( profile(lmer(..)) )")

## various scale options
xyplot(tpr.fm1,scale=list(x=list(relation="same")))  ## stupid
xyplot(tpr.fm1,scale=list(y=list(relation="same")))
xyplot(tpr.fm1,scale=list(y=list(relation="same"),tck=0))

##
expect_error(xyplot(tpr.fm1,conf=50),"must be strictly between 0 and 1")

### end {profile_plots.R}

fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)

## 0.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
##
system.time( tpr <- profile(fm01ML) )

## test all combinations of 'which', including plots (but don't show plots)
wlist <- list(1:3,1:2,1,2:3,2,3,c(1,3))
invisible(lapply(wlist,function(w) xyplot(profile(fm01ML,which=w))))

(confint(tpr) -> CIpr)
print(xyplot(tpr))
##  comparing against lme4a reference values -- but lme4 returns sigma
## rather than log(sigma)
stopifnot(dim(CIpr) == c(3,2),
          all.equal(unname(CIpr[".sigma",]),exp(c(3.64362, 4.21446)), tolerance=1e-6),
          all.equal(unname(CIpr["(Intercept)",]),c(1486.451500,1568.548494)))

options(oo)# warnings allowed ..

## fixed-effect profiling with vector RE
data(Pastes)
fmoB <- lmer(strength ~ 1 + (cask | batch), data=Pastes,
             control = lmerControl(optimizer = "bobyqa"))
(pfmoB <- profile(fmoB, which = "beta_", alphamax=.001))
xyplot(pfmoB)# nice and easy ..

summary(
    fm <- lmer(strength ~ 1 + (cask | batch), data=Pastes,
               control = lmerControl(optimizer = "nloptwrap",
                                     calc.derivs= FALSE))
)

ls.str(environment(nloptwrap))# showing *its* defaults

pfm <- profile(fm, which = "beta_", alphamax=.001) # 197 warnings for "nloptwrap"
summary(warnings())
str(pfm) # only 3 rows, .zeta = c(0, NaN, Inf) !!!
try( xyplot(pfm) ) ## FIXME or rather the profiling or rather the "wrap on nloptr"

(testLevel <- lme4:::testLevel())
if(testLevel > 2) {

    ## 2D profiles
    fm2ML <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin, REML=0)
    system.time(pr2 <- profile(fm2ML)) # 5.2 sec, 2018-05: 2.1"
    (confint(pr2) -> CIpr2)

    lme4a_CIpr2 <-
        structure(c(0.633565787613112, 1.09578224011285, -0.721864513060904,
                    21.2666273835452, 1.1821039843372, 3.55631937954106, -0.462903300019305,
                    24.6778176174587), .Dim = c(4L, 2L), .Dimnames = list(c(".sig01",
                           ".sig02", ".lsig", "(Intercept)"), c("2.5 %", "97.5 %")))
    lme4a_CIpr2[".lsig",] <- exp(lme4a_CIpr2[".lsig",])

    stopifnot(all.equal(unname(CIpr2),unname(lme4a_CIpr2),tolerance=1e-6))

    print(xyplot(pr2, absVal=0, aspect=1.3, layout=c(4,1)))
    print(splom(pr2))

    gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                 data = cbpp, family = binomial)

    ## GLMM profiles
    system.time(pr4 <- profile(gm1))  ## ~ 10 seconds

    pr4.3 <- profile(gm1,which=3)
    xyplot(pr4,layout=c(5,1),as.table=TRUE)

    splom(pr4) ## used to fail because of NAs

    nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
                 Orange, start = c(Asym = 200, xmid = 725, scal = 350))
    if (FALSE) {
        ## not working yet: detecting (slightly) lower deviance; not converging in 10k
        pr5 <- profile(nm1,which=1,verbose=1,maxmult=1.2)
        xyplot(.zeta~.focal|.par,type=c("l","p"),data=lme4:::as.data.frame.thpr(pr5),
               scale=list(x=list(relation="free")),
               as.table=TRUE)
    }
}  ## testLevel > 2

if (testLevel > 3) {
    fm3ML <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=FALSE)
    ## ~ 4 theta-variables (+ 2 fixed), 19 seconds | 2018-05: 7.4"
    print(system.time(pr3 <- profile(fm3ML)))
    print(xyplot(pr3))
    print(splom(pr3))

    if (testLevel > 4) {
      if(requireNamespace("mlmRev")) {
        data("Contraception", package="mlmRev")
        ## fit already takes ~ 3 sec (2018-05)
        fm2 <- glmer(use ~ urban+age+livch + (urban|district), Contraception, binomial)
        print(system.time(pr5 <- profile(fm2,verbose=10))) # 2018-05: 462 sec = 7'42"
        ## -> 5 warnings notably "non-monotonic profile for .sig02" (the RE's corr.)
        print(xyplot(pr5))
      }
    }  ## testLevel > 4

}  ## testLevel > 3

library("parallel")
if (detectCores()>1) {

    p0 <- profile(fm1, which="theta_")
    ## http://stackoverflow.com/questions/12983137/how-do-detect-if-travis-ci-or-not
    travis <- nchar(Sys.getenv("TRAVIS")) > 0
    if(.Platform$OS.type != "windows" && !travis) {
        prof01P <- profile(fm1, which="theta_", parallel="multicore", ncpus=2)
        stopifnot(all.equal(p0,prof01P))
    }

    ## works in Solaris from an interactive console but not ???
    ##   via R CMD BATCH

    if (Sys.info()["sysname"] != "SunOS" && !travis) {
        prof01P.snow <- profile(fm1, which="theta_", parallel="snow", ncpus=2)
        stopifnot(all.equal(p0,prof01P.snow))
    }
}

## test profile/update from within functions
foo <- function() {
    gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                 data = cbpp, family = binomial)
    ## return
    profile(gm1, which="theta_")
}
stopifnot(inherits(foo(), "thpr"))
} ## testLevel>1