source("ADMButils.s")  # Load utility functions
source("nmix_ADMB_funs.R")  # Load utility functions
source("../DATA/simNmix.R") # Simulate dataset

nID = apply(gID,1,function(x) length(unique(x)))
ID = apply(gID,1,unique)

IDind = t(apply(gID,1,find_index))

dat = list(R=R, T=T, S=21, nG=nG, N=0:20, nID=nID,ID=ID,IDind=IDind, x=x, y=y)
dat_write("nmix",dat)
pin_write("nmix",list(log_lambda=log(lambda),p0=0,p1=0,log_sigma=0,u=rep(0,nG)))

system("admb -r nmix")
st1 = system.time(system("./nmix -shess"))	# Run ADMB program

std = read.table("nmix.std",skip=1,header=F)
par = par_read("nmix.par")

coefs = std[,3]
names(coefs) = as.character(std[,2])
sd = std[,4]
names(sd) = as.character(std[,2])
confint = cbind(coefs-2*sd,coefs-2*sd)
colnames(confint) = c("2.5 %","97.5 %")

results <- list(obj=-par$loglik,  ## objective function/Neg log likelihood (NULL for MCMC)
                coef=coefs,  ## coefficients
                sd=sd,
                confint=confint,
                time=st1["elapsed"],
                terminfo=list(maxgrad=par$gradient) ## MLE termination info
                )

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

system("rm -f nmix.b01 nmix.bgs nmix.cpp nmix.eva nmix.log admodel.cov nmix.shess nmix.bar nmix.htp nmix.luu nmix.p01 nmix.rhes")
system("rm -f admodel.dep eigv.rpt  fmin.log  hessian.bin variance admodel.hes hesscheck nmix.cor  nmix.o nmix.exe")