library(aster)
library(trust)

data(radish)

pred <- c(0,1,2)
fam <- c(1,3,2)

### need object of type aster to supply to penmlogl and pickle

aout <- aster(resp ~ varb + fit : (Site * Region + Block + Pop),
    pred, fam, varb, id, root, data = radish)

### model matrices for fixed and random effects

modmat.fix <- model.matrix(resp ~ varb + fit : (Site * Region),
    data = radish)
modmat.blk <- model.matrix(resp ~ 0 + fit:Block, data = radish)
modmat.pop <- model.matrix(resp ~ 0 + fit:Pop, data = radish)

rownames(modmat.fix) <- NULL
rownames(modmat.blk) <- NULL
rownames(modmat.pop) <- NULL

idrop <- match(aout$dropped, colnames(modmat.fix))
idrop <- idrop[! is.na(idrop)]
modmat.fix <- modmat.fix[ , - idrop]

nfix <- ncol(modmat.fix)
nblk <- ncol(modmat.blk)
npop <- ncol(modmat.pop)

### try penmlogl

theta.start <- c(1, 1)

alpha.start <- aout$coefficients[match(colnames(modmat.fix),
    names(aout$coefficients))]
parm.start <- c(alpha.start, rep(0, nblk + npop))

tout <- trust(objfun = penmlogl, parm.start, rinit = 1, rmax = 10,
    sigma = theta.start, fixed = modmat.fix,
    random = list(modmat.blk, modmat.pop), obj = aout)
tout$converged

### crude estimate of variance components

eff.blk <- tout$argument[seq(nfix + 1, nfix + nblk)]
eff.pop <- tout$argument[seq(nfix + nblk + 1, nfix + nblk + npop)]
theta.crude <- sqrt(c(var(eff.blk), var(eff.pop)))

### now for the testing

pout <- penmlogl(tout$argument, theta.crude, fixed = modmat.fix,
    random = list(modmat.blk, modmat.pop), obj = aout)
names(pout)

eff.fix <- tout$argument[seq(1, nfix)]
eff.blk <- tout$argument[seq(nfix + 1, nfix + nblk)]
eff.pop <- tout$argument[seq(nfix + nblk + 1, nfix + nblk + npop)]
eff.blk <- theta.crude[1] * eff.blk
eff.pop <- theta.crude[2] * eff.pop

beta <- c(eff.fix, eff.blk, eff.pop)
scalevec <- rep(c(1, theta.crude), times = c(nfix, nblk, npop))
penaltyvec <- rep(c(0, 1, 1), times = c(nfix, nblk, npop))

modmat <- cbind(modmat.fix, modmat.blk, modmat.pop)
mout <- mlogl(beta, pred, fam, aout$x, aout$root, modmat, deriv = 2)

all.equal(mout$hessian * outer(scalevec, scalevec), pout$mlogl.hessian)
all.equal(tout$argument, pout$argument)
all.equal(scalevec, pout$scale)

all.equal(pout$value,
    mout$value + sum(penaltyvec * pout$argument^2) / 2)
all.equal(pout$gradient,
    mout$gradient * scalevec + penaltyvec * pout$argument)
all.equal(pout$hessian, pout$mlogl.hessian + diag(penaltyvec))

epsilon <- 1e-7
mygradient <- 0 * pout$gradient
myhessian <- 0 * pout$hessian
for (i in seq(along = mygradient)) {
    arg <- tout$argument
    arg[i] <- arg[i] + epsilon
    pout.epsilon <- penmlogl(arg, theta.crude, fixed = modmat.fix,
        random = list(modmat.blk, modmat.pop), obj = aout)
    mygradient[i] <- (pout.epsilon$value - pout$value) / epsilon
    myhessian[i, ] <- (pout.epsilon$gradient - pout$gradient) / epsilon
}
all.equal(mygradient, pout$gradient, tol = 100 * epsilon)
all.equal(myhessian, pout$hessian, tol = 100 * epsilon)

### penmlogl2

foo <- pout$argument
alpha2 <- foo[seq(1, nfix)]
parm2 <- foo[- seq(1, nfix)]
pout2 <- penmlogl2(parm2, alpha2, theta.crude, fixed = modmat.fix,
    random = list(modmat.blk, modmat.pop), obj = aout)
identical(pout$value, pout2$value)
identical(pout$argument, pout2$argument)
identical(pout$scale[- seq(1, nfix)], pout2$scale)
identical(pout$mlogl.hessian, pout2$mlogl.hessian)
identical(pout$mlogl.gradient, pout2$mlogl.gradient)
idx <- seq(nfix + 1, nfix + nblk + npop)
identical(pout2$gradient, pout$gradient[idx])
foom <- pout$hessian
foom <- foom[idx, ]
foom <- foom[ , idx]
identical(pout2$hessian, foom)