library(R2admb)
setup_admb()
## Read in the data
load("../DATA/Owls.rda")
sessionInfo()
## Organize the data
## these need to be defined in advance in order to be checked properly ...

## NB reported standard errors for nest effect are for sds, not variances ...
LBrood <- log(Owls$BroodSize)
Owls$ArrivalTime <- scale(Owls$ArrivalTime,center=TRUE,scale=FALSE)
mmf <- model.matrix(~(FoodTreatment+ArrivalTime)*SexParent,
                    data=Owls)
mmr <- model.matrix(~Nest-1, data=Owls)
response <- Owls$SiblingNegotiation
nf <- ncol(mmf)
nobs <- nrow(Owls)
nnests <- length(levels(Owls$Nest))
regdata=list(nobs=nrow(Owls),
     nnests=length(levels(Owls$Nest)),
     nf=ncol(mmf),
     fdesign=mmf,
     rdesign=mmr,
     obs=response,
     Lbrood=LBrood)
write_dat("owls", L=regdata)

regparams=list(fcoeffs=rep(0,ncol(mmf)),
     logit_p=0,
     sigma=1,
	rcoeffs=rep(0.001,length(levels(Owls$Nest))))
write_pin("owls", L=regparams)

compile_admb("owls", re=TRUE)

st1 <- system.time(run_admb("owls",extra.args="-noinit -nox -l1 40000000 -nl1 40000000 -ilmn 5"))["elapsed"]
#st1 <- system.time(run_admb("owls"))["elapsed"]

st1 <- unname(st1)
fit_admb <- read_admb("owls")
names(fit_admb$coeflist$fcoeffs) <-
  names(fit_admb$coefficients)[1:ncol(mmf)] <- colnames(mmf)

ev <- try(eigen(solve(vcov(fit_admb)))$value)
eratio <- if (inherits(ev,"try-error")) NA else min(ev)/max(ev)

c1 <- coef(fit_admb)[1:8]
c2 <- confint(fit_admb,method="quad")[1:8,]
## rescale from sd to var
c1[8] <- c1[8]^2
c2[8] <- c2[8]^2
regular <- list(obj=-logLik(fit_admb),  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=c1,  ## coefficients
                sd=stdEr(fit_admb)[1:8],
                confint.quad=c2,
                time=c(fit=st1),
                terminfo=c(maxgrad=fit_admb$maxgrad,eratio=eratio) ## MLE termination info
                )
clean_admb("owls")

##Now with a separable function
sepdat=list(nobs=nrow(Owls),
     nnests=length(levels(Owls$Nest)),
     nf=ncol(mmf),
     fdesign=mmf,
     nests=as.integer(Owls$Nest),
     obs=response,
     Lbrood=LBrood)
write_dat("owls_sep", L=sepdat)

seppars=list(fcoeffs=rep(0,ncol(mmf)),
     logit_p=0,
     sigma=1,
	rcoeffs=rep(0.001,length(levels(Owls$Nest))))
write_pin("owls_sep", L=seppars)
compile_admb("owls_sep", safe=FALSE, re=TRUE, verbose=FALSE)
st1_sep=system.time(run_admb("owls_sep",verbose=FALSE,extra.args="-shess -noinit -nox -l1 40000000 -nl1 40000000 -ilmn 5"))["elapsed"]
st1_sep <- unname(st1_sep)


fit_admb_sep <- read_admb("owls_sep")

c1 <- coef(fit_admb_sep)[1:8]
c2 <- confint(fit_admb_sep,method="quad")[1:8,]
## rescale from sd to var
c1[8] <- c1[8]^2
c2[8] <- c2[8]^2

ev <- try(eigen(solve(vcov(fit_admb_sep)))$value)
eratio <- if (inherits(ev,"try-error")) NA else min(ev)/max(ev)

separable <- list(obj=-logLik(fit_admb_sep),  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=c1,  ## coefficients
                sd=stdEr(fit_admb_sep)[1:8],
                confint.quad=c2,
                time=c(fit=st1_sep),
                terminfo=c(maxgrad=fit_admb_sep$maxgrad,eratio=eratio) ## MLE termination info
                )
clean_admb("owls_sep")
results=list(regular=regular, separable=separable)

save("fit_admb","fit_admb_sep","st1","st1_sep",
     file="fit_fits.RData")
save("results",file="fit.RData")