#########################################################################
##       R PROGRAM: PSE_estimate_BMA-6-a.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: Estimate population sizes from simulated samples
##                  using Bayesian loglinear model averaging  as
##                  implemented in the dga package
##
##                  Samples:
##                  Mht.10K.200, Mht.1K.025, Mht.1K.050, Mht.1K.100,
##                  Mht.1K.150, Mht.1K.200
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
##
##            DATE: 2021-02-26
##
##      DISCLAIMER: Although this program has been used by the Centers
##                  for Disease Control & Prevention (CDC), no warranty,
##                  expressed or implied, is made by the CDC or the U.S.
##                  Government as to the accuracy and functioning of the
##                  program and related program material nor shall the
##                  fact of distribution constitute any such warranty,
##                  and no responsibility is assumed by the CDC in
##                  connection therewith.
##
#########################################################################
library(HDInterval)
library(tidyr)
library(dga)
data(graphs3)
data(graphs4)
data(graphs5)
#########################################################################
## 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")
attach(dataname)
#########################################################################
## BMA loglinear model fitting to independent Mht capture
## probabilities
#########################################################################
## (t.start <- Sys.time())
## test <- bma.estimate(data = Mht.1K.200, events = 3:5,
##                      priorNmax = 10000, nreps = 1)
## Sys.time() - t.start
## print(test)
## Runs:
t.start0 <- t.start <- Sys.time()
bma.Mht.10K.200 <- bma.estimate(data = Mht.10K.200, events = 3:5,
                               priorNmax = 100000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
bma.Mht.1K.025 <- bma.estimate(data = Mht.1K.025, events = 3:5,
                               priorNmax = 10000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
bma.Mht.1K.050 <- bma.estimate(data = Mht.1K.050, events = 3:5,
                               priorNmax = 10000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
bma.Mht.1K.100 <- bma.estimate(data = Mht.1K.100, events = 3:5,
                               priorNmax = 10000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time())
bma.Mht.1K.150 <- bma.estimate(data = Mht.1K.150, events = 3:5,
                               priorNmax = 10000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
(t.start <- Sys.time() )
bma.Mht.1K.200 <- bma.estimate(data = Mht.1K.200, events = 3:5,
                               priorNmax = 10000, nreps = 400,
                               progress = FALSE)
Sys.time() - t.start
## Total elapsed time:
Sys.time() - t.start0

#########################################################################
## Save the results
#########################################################################
save(bma.Mht.10K.200, bma.Mht.1K.025, bma.Mht.1K.050, bma.Mht.1K.100,
     bma.Mht.1K.150, bma.Mht.1K.200,
     file = file.path(datapath, "BMA_Mht_estimates-a.Rdata"))
cgwtools::lsdata(file.path(datapath, "BMA_Mht_estimates-a.Rdata"))
################################   END of FILE   ########################
