library(R2admb)  ## requires version >= 0.7.2
library(coda)
source("tadpole_ADMB_funs.R")

sessionInfo()

setup_admb()
ReedfrogSizepred <- read.table("../DATA/tadpole.dat",header=TRUE)

tadpole_ADMB_setup()
st1 <- system.time(tfit_admb <- tadpole_ADMB_fit())["elapsed"]

rdat <- c(list(nobs=nrow(ReedfrogSizepred),
               nexposed=rep(10,nrow(ReedfrogSizepred))),
          ReedfrogSizepred)
rstart <- list(c=0.45,d=13,g=1)
rtype <- c(TBL="numeric",Kill="numeric")
rbounds <- list(c=c(0,1),d=c(0,50),g=c(-1,25))

st2 <- system.time(tfit_admb_prof <- do_admb("tadpole0",
                                             data=rdat,
                                             data_type=rtype,
                                             params=rstart,
                                             bounds=rbounds,
                                             run.opts=run.control(checkparam="write",
                                               checkdata="write",clean=FALSE),
                                             profile=TRUE,
                                             admb_errors="warn",
                                             profpars=c("c","d","g")))

st3 <- system.time(tfit_admb_mcmc <- do_admb("tadpole0",
                                             data=rdat,
                                             data_type=rtype,
                                             params=rstart,
                                             bounds=rbounds,
                                             run.opts=run.control(checkparam="write",
                                               checkdata="write",clean=TRUE),
                                             mcmc=TRUE,
                                             mcmc.opts=mcmc.control(mcmcpars=c("c","d","g"))))

save(tfit_admb,tfit_admb_prof,tfit_admb_mcmc,file="tadpole_ADMB_fits.RData")
## ugh: re-inverting vcov to get hessian
ev <- eigen(solve(vcov(tfit_admb)))$value
eratio <- min(ev)/max(ev)

  
results <- list(obj=-logLik(tfit_admb),  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=coef(tfit_admb),  ## coefficients
                sd=stdEr(tfit_admb),
                confint.quad=confint(tfit_admb,method="quad"),
                confint.profile=confint(tfit_admb_prof,method="profile"),
                confint.mcmc=confint(tfit_admb_mcmc,method="quantile"),
                time=c(fit=st1,
                  profile=st2,
                  mcmc=st3),  
                terminfo=c(maxgrad=tfit_admb$maxgrad,eratio=eratio), ## MLE termination info
                convinfo=list(effsize=effectiveSize(tfit_admb_mcmc$mcmc),
                  geweke=geweke.diag(tfit_admb_mcmc$mcmc))## MCMC convergence info
                )


save("results",file="fitx.RData")