#########################################################################
##       R PROGRAM: PSE_estimate_LLM-2-a.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: Estimate population sizes from simulated samples
##                  using loglinear models as implemented in the
##                  Rcapture package
##
##                  Mt, Mbh and Mht samples
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
##
##            D1ATE: 2021-03-02
#########################################################################
library(Rcapture)
#########################################################################
## Log session information
#########################################################################
sessionInfo()
#########################################################################
## Define file paths
#########################################################################
basepath <- file.path(Sys.getenv("PROJ"), "PSE/PSEsim")
workpath <- file.path(basepath, "R")
datapath <- file.path(basepath, "data")
setwd(workpath)
#########################################################################
## Source functions and get data
#########################################################################
source(file.path(workpath, "PSE_sim_functions.R"))
dataname <- file.path(datapath, "PSEsimSamples.Rdata")
load(dataname)
#########################################################################
## loglinear model fitting to independent homogeneous capture
## probabilities
#########################################################################
## test <- llm.estimate(data = M0.1K.200, events = 2:5,
##                      nreps = 2)
## print(test)
## Mt data
t.start0 <- t.start <- Sys.time()
###
llm.Mt.10K.200 <- llm.estimate(data = Mt.10K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.1K.025 <- llm.estimate(data = Mt.1K.025, events = 2:5,
                              nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.1K.050 <- llm.estimate(data = Mt.1K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.1K.100 <- llm.estimate(data = Mt.1K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.1K.150 <- llm.estimate(data = Mt.1K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.1K.200 <- llm.estimate(data = Mt.1K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
###
(t.start <- Sys.time() )
llm.Mt.20K.025 <- llm.estimate(data = Mt.20K.025, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.20K.050 <- llm.estimate(data = Mt.20K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.20K.100 <- llm.estimate(data = Mt.20K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.20K.150 <- llm.estimate(data = Mt.20K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mt.20K.200 <- llm.estimate(data = Mt.20K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
## Mbh data
(t.start <- Sys.time() )
llm.Mbh.10K.200 <- llm.estimate(data = Mbh.10K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
###
(t.start <- Sys.time() )
llm.Mbh.1K.025 <- llm.estimate(data = Mbh.1K.025, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.1K.050 <- llm.estimate(data = Mbh.1K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.1K.100 <- llm.estimate(data = Mbh.1K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.1K.150 <- llm.estimate(data = Mbh.1K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.1K.200 <- llm.estimate(data = Mbh.1K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
###
(t.start <- Sys.time() )
llm.Mbh.20K.025 <- llm.estimate(data = Mbh.20K.025, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.20K.050 <- llm.estimate(data = Mbh.20K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.20K.100 <- llm.estimate(data = Mbh.20K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.20K.150 <- llm.estimate(data = Mbh.20K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mbh.20K.200 <- llm.estimate(data = Mbh.20K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
## Mht datat.start <- Sys.time()
(t.start <- Sys.time() )
llm.Mht.10K.200 <- llm.estimate(data = Mht.10K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
###
(t.start <- Sys.time() )
llm.Mht.1K.025 <- llm.estimate(data = Mht.1K.025, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.1K.050 <- llm.estimate(data = Mht.1K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.1K.100 <- llm.estimate(data = Mht.1K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.1K.150 <- llm.estimate(data = Mht.1K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.1K.200 <- llm.estimate(data = Mht.1K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
###
(t.start <- Sys.time() )
llm.Mht.20K.025 <- llm.estimate(data = Mht.20K.025, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.20K.050 <- llm.estimate(data = Mht.20K.050, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.20K.100 <- llm.estimate(data = Mht.20K.100, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.20K.150 <- llm.estimate(data = Mht.20K.150, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
llm.Mht.20K.200 <- llm.estimate(data = Mht.20K.200, events = 2:5,
                               nreps = 400, progress = FALSE)
Sys.time() - t.start
## Total elapsed time
Sys.time() - t.start0
rm(t.start, t.start0)
#########################################################################
## Save the results
#########################################################################
save(llm.Mt.10K.200,
     llm.Mt.1K.025, llm.Mt.1K.050, llm.Mt.1K.100, llm.Mt.1K.150,
     llm.Mt.1K.200,
     llm.Mt.20K.025, llm.Mt.20K.050, llm.Mt.20K.100, llm.Mt.20K.150,
     llm.Mt.20K.200,
     llm.Mbh.10K.200,
     llm.Mbh.1K.025, llm.Mbh.1K.050, llm.Mbh.1K.100, llm.Mbh.1K.150,
     llm.Mbh.1K.200,
     llm.Mbh.20K.025, llm.Mbh.20K.050, llm.Mbh.20K.100, llm.Mbh.20K.150,
     llm.Mbh.20K.200,
     llm.Mht.10K.200,
     llm.Mht.1K.025, llm.Mht.1K.050, llm.Mht.1K.100, llm.Mht.1K.150,
     llm.Mht.1K.200,
     llm.Mht.20K.025, llm.Mht.20K.050, llm.Mht.20K.100, llm.Mht.20K.150,
     llm.Mht.20K.200,
     file = file.path(datapath, "LLM_estimates-2-a.Rdata"))
cgwtools::lsdata(file.path(datapath, "LLM_estimates-2-a.Rdata"))
################################   END of FILE   ########################
