testing.R

pacman::p_load(devtools)
load_all()
#########################################################################
sample.prior.lhs <- function(n) {
  # n: the number of samples desired
  draws0 <- randomLHS(n=n,k=8)
  draws  <- data.frame( mu_e  = qlnorm(draws0[,1],log(0.05)-1/2*0.5^2,0.5),
                        mu_l  = qlnorm(draws0[,2],log(0.25)-1/2*0.5^2,0.5),
                        mu_t  = qlnorm(draws0[,3],log(0.025)-1/2*0.5^2,0.5),
                        p     = qlnorm(draws0[,4],log(0.1)-1/2*0.5^2,0.5),
                        r_l   = qlnorm(draws0[,5],log(0.5)-1/2*0.5^2,0.5),
                        rho   = qlnorm(draws0[,6],log(0.5)-1/2*0.5^2,0.5),
                        b     = qbeta(draws0[,7],2,8),
                        c     = qlnorm(draws0[,8],log(1000)-1/2*0.2^2,0.2)
  )
  return(as.matrix(draws))
}

### A LINEAR REGRESSION EXAMPLE ####
## Define a Bayesian linear regression model
li_reg<-function(pars,data)
{
  a<-pars[1] #intercept
  b<-pars[2] #slope
  sd_e<-pars[3] #error (residuals)
  if(sd_e<=0){return(NaN)}
  pred <- a + b * data[,1]
  LLK<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
  prior<- prior_reg(pars)
  return(LLK + prior)
}
## Define the Prior distributions
prior_reg<-function(pars)
{
  a<-pars[1] #intercept
  b<-pars[2] #slope
  epsilon<-pars[3] #error
  prior_a<-dnorm(a,0,100,log=TRUE) ## non-informative (flat) priors on all
  prior_b<-dnorm(b,0,100,log=TRUE) ## parameters.
  prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)
  return(prior_a + prior_b + prior_epsilon)
}

# simulate data
x<-runif(30,5,15)
y<-x+rnorm(30,0,5)
d<-cbind(x,y)
mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
                       par_names=c('a','b','epsilon'),data=d)
## For best results, run again with the previously
## adapted variance-covariance matrix.
mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
                       prop_sigma=mcmc_r$prop_sigma,par_names=c('a','b','epsilon'),data=d)
mcmc_r<-mcmc_thin(mcmc_r)
plotMH(mcmc_r)


tst = HID_markov(.v_params = name_HID_params(rep(0.5, 9)))

tst2 = HID_markov(.v_params = name_HID_params(rep(0.5, 9)), project_future = T)


v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))

tst = sample_prior_LHS(.l_params = list(v_params_names = v_params_names,
                                        v_params_dists = v_params_dists, args = args), 10)

tst2 = sample_prior_FGS(.l_params = list(v_params_names = v_params_names,
                                         v_params_dists = v_params_dists, args = args),.n_samples = 10)
tst2

tst3 = sample_prior_RGS(.l_params = list(v_params_names = v_params_names,
                                         v_params_dists = v_params_dists, args = args),.n_samples = 10)
tst3
#########################################################################
# Number of initial starting points - NM:
n_init <- 100
# Number of random samples:
n_samples <- 10
# Names and number of input parameters to be calibrated:
v_params_names <- c("p_Mets", "p_DieMets")
n_params <- length(v_params_names)
v_params_init <- samples[1:10,]
n_init <- nrow(v_params_init)
## Run Gradient-based for each starting point:
res_llk = list()
res_sse = list()
m_calib_res_llk <- m_calib_res_sse <-
  matrix(nrow = n_init, ncol = n_params + 1)
colnames(m_calib_res_llk) <- colnames(m_calib_res_sse) <-
  c(v_params_names, "Overall_fit")

for (j in 1:n_init) {
  fit_sa <- optim(par = v_params_init[j, ],
                  fn = LLK, # GOF is log likelihood
                  method = "SANN",
                  control = list(  fnscale = -1,
                                   temp = 10,
                                   tmax = 10,
                                   maxit = 1000), # maximum iterations
                  hessian = TRUE,
                  .func = CRS_markov, # model to be optimised
                  .args = NULL, # arguments to be passed to the model
                  .l_targets = l_targets, # targets passed to .gof
                  .maximise = TRUE, # .gof should maximise
                  .optim = TRUE)
  m_calib_res_llk[j, ] <- c(fit_sa$par, fit_sa$value)
}

#########################################################################
load(file.path(here::here(), "data", "CRS_targets.rda"))
v_targets_names <- c("Surv")
v_targets_dists <- c('norm')
data("CRS_targets")
Surv <- CRS_targets$Surv
l_targets <- list(
  'v_targets_names' = v_targets_names,
  'Surv' = Surv,
  'v_targets_dists' = v_targets_dists)

testing <- LLK(.func = CRS_markov,
                          .samples = tst,
                          .l_targets = l_targets)
###########################
data("CRS_targets")
Surv <- CRS_targets$Surv
v_targets_names <- c("Surv", "Surv")
v_targets_dists <- c('norm', 'norm')
v_targets_weights <- c(0.2, 0.8)
l_targets <- list('v_targets_names' = v_targets_names, 'Surv' = Surv, 'v_targets_dists' = v_targets_dists, 'v_targets_weights' = v_targets_weights)
v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))

samples <- sample_prior_LHS(
  .l_params = list(v_params_names = v_params_names,                             v_params_dists = v_params_dists, args = args), .n_samples = 10000)

GOF_llik <- LLK(.func = CRS_markov, .samples = samples,
                           .l_targets = l_targets, .sample_method = "LHS")
#LLK##########################
data("CRS_targets")
Surv <- CRS_targets$Surv
v_targets_names <- c("Surv", "Surv")
v_targets_weights <- c(0.5, 0.5)
v_targets_dists <- c('norm')
# v_targets_names <- c("Surv")
# v_targets_weights <- c(1)
l_targets <-
  list('v_targets_names' = v_targets_names,
       'Surv' = Surv, 'v_targets_dists' = v_targets_dists,
       'v_targets_weights' = v_targets_weights)
v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))
l_params <- list(v_params_names = v_params_names,                             v_params_dists = v_params_dists, args = args)

samples <- sample_prior_LHS(
  .l_params = l_params, .n_samples = 10000)

GOF_llik1 <- LLK(.func = CRS_markov, .samples = samples,
                            .l_targets = l_targets,
                            .sample_method = "LHS")
# GOF_llik2 <- LLK(.func = CRS_markov, .samples = samples,
#                             .l_targets = l_targets, .optim = TRUE)
# GOF_llik4 <- LLK(.func = CRS_markov, .samples = samples,
#                             .l_targets = l_targets, .optim = TRUE)
#wSSE_GOF##########################
data("CRS_targets")
Surv <- CRS_targets$Surv
v_targets_names <- c("Surv", "Surv")
v_targets_weights <- c(0.5, 0.5)
v_targets_dists <- c('norm')
# v_targets_names <- c("Surv")
# v_targets_weights <- c(1)
l_targets <-
  list('v_targets_names' = v_targets_names,
       'v_targets_weights' = v_targets_weights,
       'Surv' = Surv)
v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))

# samples <- sample_prior_LHS(
#   .l_params = list(v_params_names = v_params_names,                             v_params_dists = v_params_dists, args = args), .n_samples = 10000)

GOF_wsse1 <- wSSE_GOF(.func = CRS_markov, .samples = samples,
                      .l_targets = l_targets, .sample_method = "LHS")
# GOF_wsse2 <- wSSE_GOF(.func = CRS_markov, .samples = samples[1:2,],
#                       .l_targets = l_targets)
# GOF_wsse5 <- wSSE_GOF(.func = CRS_markov, .samples = samples,
#                       .l_targets = l_targets, .optim = TRUE)
# GOF_wsse4 <- wSSE_GOF(.func = CRS_markov, .samples = samples,
#                       .l_targets = l_targets)
# GOF_wsse5 <- wSSE_GOF(.func = CRS_markov, .samples = samples,
#                       .l_targets = l_targets)
compare(GOF_wsse, GOF_wsse2)
###########################

tsts = function(...) {
  dots = list(...)
  # moz = NULL
  moz = dots[['moz']]
  if(is.null(moz))
    cat('I can find the default')
  cat(moz)
}
tsts(2, 3)
`if`(0, "test", "not")
##########################################################
DEoptim::DEoptim(
  fn = wSSE_GOF,
  lower = lb,
  upper = ub,
  control = DEoptim::DEoptim.control( # control parameters
    trace = FALSE), # printing a trace
  .func = CRS_markov, # model to be optimised
  .args = NULL, # arguments to be passed to the model
  .l_targets = l_targets, # targets passed to .gof
  .maximise = FALSE, # .gof should minimise
  .optim = TRUE)
#Directed_search####################################################
load_all()

data("CRS_targets")
Surv <- CRS_targets$Surv
v_targets_names <- c("Surv", "Surv")
v_targets_weights <- c(0.5, 0.5)
v_targets_dists <- c("norm", "norm")
# v_targets_names <- c("Surv")
# v_targets_weights <- c(1)
l_targets <-
  list('v_targets_names' = v_targets_names,
       'Surv' = Surv,
       'v_targets_dists' = v_targets_dists,
       'v_targets_weights' = v_targets_weights)
v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))
l_params <- list('v_params_names' = v_params_names,
                 'v_params_dists' = v_params_dists,
                 'args' = args,
                 'Xargs' = args)
rm(v_params_names, v_params_dists, v_targets_dists, v_targets_weights,
   v_targets_names, args)

set.seed(1)
samples <- sample_prior_LHS(.l_params = l_params,
                            .n_samples = 50)

NM_optimise_wSSE <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000)

GB_optimise_wSSE <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000)

SA_optimise_wSSE <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = list(NULL),
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

GA_optimise_wSSE <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = list(NULL),
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

NM_optimise_lLLK <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000)

GB_optimise_lLLK <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000)

SA_optimise_lLLK <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

GA_optimise_lLLK <- calibrateModel_directed(
  .l_params = l_params,
  .func = CRS_markov,
  .args = list(NULL),
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

#PSA_values##################################################

l_optim_lists <- list(GA_optimise_lLLK, GA_optimise_wSSE, GB_optimise_lLLK,
                      GB_optimise_wSSE, NM_optimise_lLLK, NM_optimise_wSSE,
                      SA_optimise_lLLK, SA_optimise_wSSE)

testing <- PSA_calib_values(.l_optim_lists = l_optim_lists,
                            .search_method = "Directed")
testing %>% transpose() %>% View()


l_optim_lists2 <- list(GOF_llik1, GOF_wsse1)

testing2 <- PSA_calib_values(.l_optim_lists = l_optim_lists2,
                             .search_method = "Random")

#Bayesian_calibration_helpers###########################################
##### prior:
log_prior(.samples = samples[1,], .l_params = l_params)
calc_log_prior <- function(.n_param = n_params, .v_params,
                           .v_params_names = v_params_names) {
  if(is.null(dim(.v_params))) { # If vector, change to matrix
    .v_params <- t(.v_params)
  }
  n_samp <- nrow(.v_params)
  colnames(.v_params) <- .v_params_names
  lprior <- rep(0, n_samp)
  for (i in 1:.n_param){
    lprior <- lprior + dunif(.v_params[, i],
                             min = lb[i],
                             max = ub[i],
                             log = TRUE)
    # ALTERNATIVE prior using beta distributions
    # lprior <- lprior + dbeta(v_params[, i],
    #                          shape1 = 1,
    #                          shape2 = 1,
    #                          log = T)
  }
  return(lprior)
}
lb <- c(p_Mets = 0.04, p_DieMets = 0.04) # lower bound
ub <- c(p_Mets = 0.16, p_DieMets = 0.12) # upper bound
v = samples %>% as.matrix()
calc_log_prior(.v_params = v, .n_param = 2,
               .v_params_names = v_params_names) ==
  log_prior(.samples = samples, .l_params = l_params)

log_prior(.samples = samples, .l_params = l_params)
calculate_prior(.samples = samples[1,], .l_params = l_params)
##### likelihood:
calc_log_lik <- function(.func = CRS_markov, .lst_targets = lst_targets,
                         .v_params, .n_target = n_target){
  if(is.null(dim(.v_params))) { # If vector, change to matrix
    .v_params <- t(.v_params)
  }
  n_samp <- nrow(.v_params)
  v_llik <- matrix(0, nrow = n_samp, ncol = .n_target)
  llik_overall <- numeric(n_samp)
  for(j in 1:n_samp) { # j=1
    jj <- tryCatch( {
      ### Run model for a given parameter set:
      model_res <- exec(.fn = .func, .v_params[j, ])

      ###  Calculate log-likelihood of model outputs to targets  ###
      # TARGET 1: Survival ("Surv")
      # log likelihood
      v_llik[j, 1] <- sum(dnorm(x = .lst_targets$Surv$value,
                                mean = model_res$Surv,
                                sd = .lst_targets$Surv$se,
                                log = TRUE))

      # TARGET 2: (if you had more...)
      # log likelihood
      # v_llik[j, 2] <- sum(dnorm(x = lst_targets$Target2$value,
      #                        mean = model_res$Target2,
      #                        sd = lst_targets$Target2$se,
      #                        log = T))

      # OVERALL
      llik_overall[j] <- sum(v_llik[j, ])
    }, error = function(e) NA)
    if(is.na(jj)) { llik_overall <- -Inf }
  } # End loop over sampled parameter sets
  # return LLIK
  return(llik_overall)
}

LLK(.samples = samples, .func = CRS_markov,
               .args = NULL, .l_targets = l_targets)
calc_log_lik(.v_params = v)
calculate_likelihood(.samples = samples, .func = CRS_markov,
                     .args = NULL, .l_targets = l_targets)
calc_likelihood(.v_params = v)
calc_likelihood(.v_params = v) ==
  calculate_likelihood(.samples = samples, .func = CRS_markov,
                       .args = NULL, .l_targets = l_targets)
##### posterior:
calc_log_post <- function(.v_params, .target = lst_targets) {
  # Call log-likelihood function:
  LLK <- calc_log_lik(.v_params = .v_params,
                                 .lst_targets = .target)
  # Call log-prior function:
  lprior <- calc_log_prior(.v_params = .v_params)
  # Compute log-posterior:
  lpost <- LLK + lprior

  return(lpost)
}
calc_log_post(.v_params = v[1,], .target = lst_targets)
exp(calc_log_post(.v_params = v[1,], .target = lst_targets))

log_posterior(.samples = samples, .func = CRS_markov, .args = NULL,
              .l_targets = l_targets, .l_params = l_params)
calculate_posterior(.samples = samples, .func = CRS_markov,
                    .args = NULL, .l_targets = l_targets,
                    .l_params = l_params)

#Bayesian_calibration##################################################
data("CRS_targets")
Surv <- CRS_targets$Surv
v_targets_names <- c("Surv", "Surv")
v_targets_weights <- c(0.5, 0.5)
v_targets_dists <- c("norm", "norm")
# v_targets_names <- c("Surv")
# v_targets_weights <- c(1)
l_targets <-
  list('v_targets_names' = v_targets_names,
       'Surv' = Surv,
       'v_targets_dists' = v_targets_dists,
       'v_targets_weights' = v_targets_weights)
v_params_names <- c("p_Mets", "p_DieMets")
v_params_dists <- c("unif", "unif")
args <- list(list(min = 0.04, max = 0.16),
             list(min = 0.04, max = 0.12))
l_params <- list('v_params_names' = v_params_names,
                 'v_params_dists' = v_params_dists,
                 'args' = args,
                 'Xargs' = args)
rm(v_params_names, v_params_dists, v_targets_dists, v_targets_weights,
   v_targets_names, args)

set.seed(1)
samples <- sample_prior_LHS(.l_params = l_params,
                            .n_samples = 1000)

test_Bayesian = calibrateModel_beyesian(
  .b_method = 'SIR', .func = CRS_markov, .args = NULL,
  .l_targets = l_targets, .l_params = l_params, .samples = samples)

test_Bayesian2 = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = CRS_markov, .args = NULL,
  .l_targets = l_targets, .l_params = l_params, .samples = samples,
  .n_resample = 1000)

#CRS_model################################################################
library(devtools)
load_all()

samples_CRS_data <- sample_prior_LHS(
  .n_samples = 1000,
  .l_params = CRS_data$l_params)
samples1_CRS_data <- sample_prior_FGS(
  .n_samples = 5,
  .l_params = CRS_data$l_params)
samples2_CRS_data <- sample_prior_RGS(
  .n_samples = 1000,
  .l_params = CRS_data$l_params)

GOF_wsse_CRS_model <- wSSE_GOF(
  .func = CRS_markov,
  .optim = FALSE,
  .args = NULL,
  .samples = samples_CRS_data,
  .l_targets = CRS_data$l_targets,
  .sample_method = "LHS")

GOF_llik_CRS_model <- LLK_GOF(
  .func = CRS_markov, .optim = FALSE,
  .args = NULL,
  .samples = samples_CRS_data,
  .l_targets = CRS_data$l_targets,
  .sample_method = "LHS")

samples <- sample_prior_LHS(.n_samples = 100,
                            .l_params = CRS_data$l_params)

NM_optimise_wSSE_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000)

GB_optimise_wSSE_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000)

SA_optimise_wSSE_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

GA_optimise_wSSE_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

NM_optimise_LLK_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000)

GB_optimise_LLK_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000)

SA_optimise_LLK_CRS_model <- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

GA_optimise_LLK_CRS_model<- calibrateModel_directed(
  .l_params = CRS_data$l_params,
  .func = CRS_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = CRS_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

samples <- sample_prior_LHS(.n_samples = 1000,
                            .l_params = CRS_data$l_params)

SIR_CRS_data = calibrateModel_beyesian(
  .b_method = 'SIR', .func = CRS_markov,
  .args = NULL,
  .l_targets = CRS_data$l_targets,
  .l_params = CRS_data$l_params, .samples = samples)

set.seed(1) # Function crashes on set.seed(1)
IMIS2_CRS_data = calibrateModel_beyesian2(
  .b_method = 'IMIS', .func = CRS_markov,
  .args = NULL,
  .l_targets = CRS_data$l_targets,
  .l_params = CRS_data$l_params,
  .n_resample = 1000)

rm(likelihood, prior, sample.prior)
set.seed(1) # Function crashes on set.seed(1)
IMIS_CRS_data = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = CRS_markov,
  .args = NULL,
  .l_targets = CRS_data$l_targets,
  .l_params = CRS_data$l_params,
  .n_resample = 1000)

l_optim_lists_CRS_data <- list(
  GA_optimise_LLK_CRS_model, GA_optimise_wSSE_CRS_model,
  GB_optimise_LLK_CRS_model, GB_optimise_wSSE_CRS_model,
  NM_optimise_LLK_CRS_model, NM_optimise_wSSE_CRS_model,
  SA_optimise_LLK_CRS_model, SA_optimise_wSSE_CRS_model)

l_optim_lists_CRS_data2 <- list(
  GOF_wsse_CRS_model, GOF_llik_CRS_model)

l_optim_lists_CRS_data3 <- list(
  SIR_CRS_data, IMIS_CRS_data, IMIS2_CRS_data)

PSA_values_CRS_markov <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_CRS_data,
  .search_method = 'Directed',
  .PSA_runs = 1000,
  .l_params = CRS_data$l_params)

PSA_values_CRS_markov2 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_CRS_data2,
  .search_method = 'Random',
  .PSA_runs = 1000,
  .l_params = CRS_data$l_params)

PSA_values_CRS_markov3 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_CRS_data3,
  .search_method = 'Bayesian',
  .PSA_runs = 1000,
  .l_params = CRS_data$l_params)

#Run_PSA#############################################################
# calib_values_ <- c(PSA_values_HID_data,
#                    PSA_values_HID_data2,
#                    PSA_values_HID_data3)
#
# PSA_results <- run_PSA(
#   .func_ = HID_markov,
#   .args_ = NULL,
#   .PSA_calib_values_ = c(PSA_values_HID_data,
#                          PSA_values_HID_data2,
#                          PSA_values_HID_data3),
#   .PSA_unCalib_values_ = NULL)

PSA_results_CRS_markov <- run_PSA(
  .func_ = CRS_markov,
  .args_ = NULL,
  .PSA_calib_values_ = c(PSA_values_CRS_markov,
                         PSA_values_CRS_markov2,
                         PSA_values_CRS_markov3),
  .PSA_unCalib_values_ = NULL)

















#HID_model#############################################################
l_prior <- function(par_vector) {
  # par_vector: a vector (or matrix) of model parameters (omits c)
  if(is.null(dim(par_vector))) par_vector <- t(par_vector)
  lprior <- rep(0,nrow(par_vector))
  lprior <- lprior+dlnorm(par_vector[,1],log(0.05 )-1/2*0.5^2,0.5,log=TRUE)    # mu_e
  lprior <- lprior+dlnorm(par_vector[,2],log(0.25 )-1/2*0.5^2,0.5,log=TRUE)    # mu_l
  lprior <- lprior+dlnorm(par_vector[,3],log(0.025)-1/2*0.5^2,0.5,log=TRUE)    # mu_t
  lprior <- lprior+dlnorm(par_vector[,4],log(0.1  )-1/2*0.5^2,0.5,log=TRUE)    # p
  lprior <- lprior+dlnorm(par_vector[,5],log(0.5  )-1/2*0.5^2,0.5,log=TRUE)    # r_l
  lprior <- lprior+dlnorm(par_vector[,6],log(0.5  )-1/2*0.5^2,0.5,log=TRUE)    # rho
  lprior <- lprior+dbeta( par_vector[,7],2,8,log=TRUE)                         # b
  return(lprior)
}

tooot = samples_HID2_data[1:5,]
tooot2 = backTransform(.t_data_ = as_tibble(samples_HID2_data[1:5,]), .l_params_ = HID_data2$l_params)

log_prior(.samples = tooot2[1,], .l_params = HID_data$l_params)
l_prior(par_vector = as.matrix(tooot2[1,]))
log_prior(.samples = tooot[1,], .l_params = HID_data2$l_params,
          .transform = TRUE)

microbenchmark::microbenchmark(
  log_prior_(.samples = samples_HID2_data, .l_params = HID_data2$l_params,
             .transform = TRUE),
  l_prior(par_vector = as.matrix(samples_HID2_data)))

mod <- function(par_vector,project_future=FALSE) {
  # par_vector: a vector of model parameters
  # project_future: TRUE/FALSE, whether to project future outcomes for policy comparison
  pop_size   <- 1e6             # population size hard-coded as 1 million
  mu_b       <- 0.015           # background mortality rate hard-coded as 0.015
  mu_e       <- par_vector[1]   # cause-specific mortality rate with early-stage disease
  mu_l       <- par_vector[2]   # cause-specific mortality rate with late-stage disease
  mu_t       <- par_vector[3]   # cause-specific mortality rate on treatment
  p          <- par_vector[4]   # transtion rate from early to late-stage disease
  r_l <- r_e <- par_vector[5]   # rate of uptake onto treatment (r_l = late-stage disease;r_e = early-stage disease)
  rho        <- par_vector[6]   # effective contact rate
  b          <- par_vector[7]   # fraction of population in at-risk group
  c          <- par_vector[8]   # annual cost of treatment

  ######## Prepare to run model ###################
  n_yrs    <- if(project_future) { 51 } else { 30 }  # no. years to simulate (30 to present, 51 for 20 year analytic horizon)
  sim      <- if(project_future) { 1:2 } else { 1 }  # which scenarios to simulate: 1 = base case, 2 = expanded treatment access
  v_mu     <- c(0,0,mu_e,mu_l,mu_t)+mu_b             # vector of mortality rates
  births   <- pop_size*mu_b*c(1-b,b)                 # calculate birth rate for equilibrium population before epidemic
  init_pop <- pop_size*c(1-b,b-0.001,0.001,0,0,0)    # creates starting vector for population
  trace    <- matrix(NA,12*n_yrs,6)                  # creates a table to store simulation trace
  colnames(trace) <- c("N","S","E","L","T","D")
  results  <- list()                                 # creates a list to store results

  ######## Run model ###################
  for(s in sim) {
    P0 <- P1 <- init_pop
    for(m in 1:(12*n_yrs)) {
      lambda    <- rho*sum(P0[3:4])/sum(P0[2:5]) # calculates force of infection
      P1[1:2]   <- P1[1:2]+births/12             # births
      P1[-6]    <- P1[-6]-P0[-6]*v_mu/12         # deaths: N, S, E, L, T, to D
      P1[6]     <- P1[6]+sum(P0[-6]*v_mu/12)     # deaths:N, S, E, L, T, to D
      P1[2]     <- P1[2]-P0[2]*lambda/12         # infection: S to E
      P1[3]     <- P1[3]+P0[2]*lambda/12         # infection: S to E
      P1[3]     <- P1[3]-P0[3]*p/12              # progression: E to L
      P1[4]     <- P1[4]+P0[3]*p/12              # progression: E to L
      P1[4]     <- P1[4]-P0[4]*r_l/12            # treatment uptake: L to T
      P1[5]     <- P1[5]+P0[4]*r_l/12            # treatment uptake: L to T
      if(s==2 & m>(12*30)) {
        P1[3]   <- P1[3]-P0[3]*r_e/12            # treatment uptake: E to T (scenario 2)
        P1[5]   <- P1[5]+P0[3]*r_e/12            # treatment uptake: E to T (scenario 2)
      }
      trace[m,] <- P0 <- P1                      # fill trace, reset pop vectors
    }
    results[[s]] <- trace                        # save results for each scenario
  }

  ######## Report results ###################
  if(project_future==FALSE) {
    ## Return calibration metrics, if project_future = FALSE
    return(list(prev = (rowSums(trace[,3:5])/rowSums(trace[,1:5]))[c(10,20,30)*12],  # Prevalence at 10,20,30 years
                surv = 1/(v_mu[3]+p)+ p/(v_mu[3]+p)*(1/v_mu[4]),                     # HIV survival without treatment
                tx   = trace[30*12,5]                                                # Treatment volume at 30 years
    ) )
  } else {
    ## Policy projections for CE analysis, if project_future = TRUE
    return(list(trace0   = results[[1]],     # Trace without expanded treatment access
                trace1   = results[[2]],     # Trace with expanded treatment access
                inc_LY   = sum(results[[2]][(30*12+1):(51*12),-6]-results[[1]][(30*12+1):(51*12),-6])/12,  # incr. LY lived with expanded tx
                inc_cost = sum(results[[2]][(30*12+1):(51*12),5]-results[[1]][(30*12+1):(51*12),5])*c/12   # incr. cost  with expanded tx
    ) )
  }
}

l_likelihood <- function(par_vector) {
  # par_vector: a vector (or matrix) of model parameters
  if(is.null(dim(par_vector))) par_vector <- t(par_vector)
  llik <- rep(0,nrow(par_vector))
  for(j in 1:nrow(par_vector)) {
    jj <- tryCatch( {
      res_j <- mod(c(as.numeric(par_vector[j,]),1))
      llik[j] <- llik[j]+sum(dbinom(c(25,75,50),500,res_j[["prev"]], log=TRUE)) # prevalence likelihood
      llik[j] <- llik[j]+dnorm(10,res_j[["surv"]],2/1.96, log=TRUE)             # survival likelihood
      llik[j] <- llik[j]+dnorm(75000,res_j[["tx"]],5000/1.96, log=TRUE)         # treatment volume likelihood
    }, error = function(e) NA)
    if(is.na(jj)) { llik[j] <- -Inf }
  }
  return(llik)
}

l_likelihood(par_vector = as.matrix(tooot2[3,]))
LLK(
  .samples = tooot2[3,],
  .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets
)

microbenchmark::microbenchmark(
  l_likelihood(par_vector = as.matrix(tooot2[3,])),
  LLK(
    .samples = tooot2[3,],
    .func = HID_markov,
    .args = NULL,
    .l_targets = HID_data$l_targets
  )
)

l_post <- function(par_vector) {
  return( l_prior(par_vector) + l_likelihood(par_vector) )
}
calculate_posterior(
  .samples = tooot2[3,],
  .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params)

l_post(as.matrix(tooot2[3,]))
log_posterior(
  .samples = tooot2[3,],
  .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params)

microbenchmark::microbenchmark(
  l_post(as.matrix(tooot2[3,])),
  log_posterior(
    .samples = tooot2[3,],
    .func = HID_markov,
    .args = NULL,
    .l_targets = HID_data$l_targets,
    .l_params = HID_data$l_params)
)

microbenchmark::microbenchmark(
  l_post(as.matrix(as_tibble(t(tooot)))),
  log_posterior(
    .samples = tooot,
    .func = HID_markov_2,
    .args = NULL,
    .l_targets = HID_data2$l_targets,
    .l_params = HID_data2$l_params)
)

# tst = HID_markov(.v_params = name_HID_params(rep(0.5, 9)))
# tst1 = HID_markov(.v_params = name_HID_params(rep(1, 9)))
# tst2 =  HID_markov()
# tst2 = HID_markov(.v_params = name_HID_params(rep(0.5, 9)), project_future = T)

# v_targets_names <- c("Prev", "Surv", "Trt_vol")
# v_targets_weights <- c(1, 1, 1)
# v_targets_dists <- c("binom", "norm", "norm")
# l_targets <-
#   list('v_targets_names' = v_targets_names,
#        'Prev' = tibble('value' = c(5/100, 15/100, 10/100), # %
#                        'se' = c(5/1000, 15/1000, 10/1000), # 10% of value
#                        'x' = c(25, 50, 75),
#                        'size' = 500,
#                        'lb' = c(3.3, 12, 7.5),
#                        'ub' = c(7.1, 18.3, 12.8)),
#        'Surv' = tibble('value' = 10,
#                        'se' = 2/1.96,
#                        'lb' = 8,
#                        'ub' = 12),
#        'Trt_vol' = tibble('value' = 75000,
#                           'se' = 5000/1.96,
#                           'lb' = 70000,
#                           'ub' = 80000),
#        'v_targets_dists' = v_targets_dists,
#        'v_targets_weights' = v_targets_weights)
# v_params_names <- c("mu_e", "mu_l", "mu_t", "p", "r_l", "rho",
#                     "b")
# v_params_dists <- c("lnorm", "lnorm", "lnorm", "lnorm", "lnorm", "lnorm",
#                     "beta")
# args <- list(list(meanlog = -3.121, sdlog = 0.5),
#              list(meanlog = -1.511, sdlog = 0.5),
#              list(meanlog = -3.814, sdlog = 0.5),
#              list(meanlog = -2.428, sdlog = 0.5),
#              list(meanlog = -0.818, sdlog = 0.5),
#              list(meanlog = -0.818, sdlog = 0.5),
#              list(shape1 = 2, shape2 = 8))
# extra_args <- list(list(min = 0.02, max = 0.12),
#                    list(min = 0.08, max = 0.59),
#                    list(min = 0.01, max = 0.06),
#                    list(min = 0.03, max = 0.24),
#                    list(min = 0.17, max = 1.18),
#                    list(min = 0.01, max = 0.06),
#                    list(min = 0.03, max = 0.48))
# l_params <- list('v_params_names' = v_params_names,
#                  'v_params_dists' = v_params_dists,
#                  'args' = args,
#                  'Xargs' = extra_args)
# rm(v_params_names, v_params_dists, v_targets_dists, v_targets_weights,
#    v_targets_names, args)
library(devtools)
load_all()

samples_HID_data <- sample_prior_LHS(.n_samples = 1000,
                                     .l_params = HID_data$l_params)
samples1_HID_data <- sample_prior_FGS(.n_samples = 5,
                                      .l_params = HID_data$l_params)
samples2_HID_data <- sample_prior_RGS(.n_samples = 50,
                                      .l_params = HID_data$l_params)

GOF_wsse_HID_data <- wSSE_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = samples_HID_data,
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")

GOF_llik_HID_data <- LLK_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = samples_HID_data,
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")

samples <- sample_prior_LHS(.n_samples = 5,
                            .l_params = HID_data$l_params)

NM_optimise_wSSE_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

GB_optimise_wSSE_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

SA_optimise_wSSE_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

GA_optimise_wSSE_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

NM_optimise_LLK_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

GB_optimise_LLK_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

SA_optimise_LLK_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

GA_optimise_LLK_HID_data <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

samples <- sample_prior_LHS(.n_samples = 1000,
                            .l_params = HID_data$l_params)

SIR_HID_data = calibrateModel_beyesian(
  .b_method = 'SIR', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params, .samples = samples)

set.seed(1) # Function crashes on set.seed(1)
IMIS2_HID_data = calibrateModel_beyesian2(
  .b_method = 'IMIS', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params,
  .n_resample = 1000)

rm(likelihood, prior, sample.prior)
set.seed(1) # Function crashes on set.seed(1)
IMIS_HID_data = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params,
  .n_resample = 1000)

l_optim_lists_HID_data <- list(
  GA_optimise_LLK_HID_data, GA_optimise_wSSE_HID_data,
  GB_optimise_LLK_HID_data, GB_optimise_wSSE_HID_data,
  NM_optimise_LLK_HID_data, NM_optimise_wSSE_HID_data,
  SA_optimise_LLK_HID_data, SA_optimise_wSSE_HID_data)

l_optim_lists_HID_data2 <- list(
  GOF_wsse_HID_data, GOF_llik_HID_data)

l_optim_lists_HID_data3 <- list(
  SIR_HID_data, IMIS_HID_data, IMIS2_HID_data)

PSA_values_HID_data <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_HID_data,
  .search_method = 'Directed',
  .PSA_runs = 1000,
  .l_params = HID_data$l_params)

PSA_values_HID_data2 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_HID_data2,
  .search_method = 'Random',
  .PSA_runs = 1000,
  .l_params = HID_data$l_params)

PSA_values_HID_data3 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_HID_data3,
  .search_method = 'Bayesian',
  .PSA_runs = 1000,
  .l_params = HID_data$l_params)

#old = options('warning.length')
#options(warning.length=8000)
#options(old)

#Rcpp microsim##############################################
# Arguments:
# v_M_1:   vector of initial states for individuals
# n.i:     number of individuals
# n.t:     total number of cycles to run the model
# v.n:     vector of health state names
# d.c:     discount rate for costs
# d.e:     discount rate for health outcome (QALYs)
# t_p:     vector containing transition probability
# u_vec:   utilities vectors
# c_vec:   costs vectors
# Trt:     are the n.i individuals receiving treatment? (scalar with a Boolean value, default is FALSE)
# seed:    starting seed number for random number generator (default is 1)
# Makes use of:
# ProbsCpp: function for the estimation of transition probabilities
# CostsRcpp2: function for the estimation of cost state values
# EffsRcpp2: function for the estimation of state specific health outcomes (QALYs)

# n.i   <- 100000                # number of simulated individuals
# n.t   <- 30                    # time horizon, 30 cycles
# v.n   <- c("H","S1","S2","D")  # the model states: Healthy (H), Sick (S1), Sicker (S2), Dead (D)
# n.s   <- length(v.n)           # the number of health states
# v.M_1 <- rep("H", n.i)         # everyone begins in the healthy state
# d.c   <- d.e <- 0.03           # equal discounting of costs and QALYs by 3%
# v.Trt <- c("No Treatment", "Treatment") # store the strategy names
#
# # Transition probabilities (per cycle)
# p.HD    <- 0.005               # probability to die when healthy
# p.HS1   <- 0.15          	     # probability to become sick when healthy
# p.S1H   <- 0.5           	     # probability to become healthy when sick
# p.S1S2  <- 0.105         	     # probability to become sicker when sick
# rr.S1   <- 3             	     # rate ratio of death in sick vs healthy
# rr.S2   <- 10            	     # rate ratio of death in sicker vs healthy
# r.HD    <- -log(1 - p.HD) 	   # rate of death in healthy
# r.S1D   <- rr.S1 * r.HD  	     # rate of death in sick
# r.S2D   <- rr.S2 * r.HD  	     # rate of death in sicker
# p.S1D   <- 1 - exp(- r.S1D)    # probability to die in sick
# p.S2D   <- 1 - exp(- r.S2D)    # probability to die in sicker
#
# # Cost and utility inputs
# c.H     <- 2000                # cost of remaining one cycle healthy
# c.S1    <- 4000                # cost of remaining one cycle sick
# c.S2    <- 15000               # cost of remaining one cycle sicker
# c.Trt   <- 12000               # cost of treatment (per cycle)
#
# u.H     <- 1                   # utility when healthy
# u.S1    <- 0.75                # utility when sick
# u.S2    <- 0.5                 # utility when sicker
# u.Trt   <- 0.95                # utility when being treated
#
# # Define starting health state, using numbers instead of characters to identify the health states:
# v_M_1 = rep(1, n.i)
# #v_M_1 = rep(c(1, 2, 3, 4), n.i/4)
#
# # Create a vector of transition probabilities:
# t_p = c(p.HD, p.HS1, p.S1H, p.S1S2, p.S1D, p.S2D)
# names(t_p) = c("p.HD", "p.HS1", "p.S1H", "p.S1S2", "p.S1D", "p.S2D")
#
# # Create a vector containing costs parameters:
# c_vec = c(c.H, c.S1, c.S2, c.Trt)
# names(c_vec) = c("c.H", "c.S1", "c.S2", "c.Trt")
#
# # Create a vector containing utilities parameters:
# u_vec = c(u.H, u.S1, u.S2, u.Trt)
# names(u_vec) = c("u.H", "u.S1", "u.S2", "u.Trt")
#
# ###
#
# ResV_no_trt_Cpp =
#   MicroSimV_Cpp(v_S_t = v_M_1, t_P = t_p, v_C = c_vec, v_U = u_vec, n_I = n.i,
#                 n_S = n.s, n_T = n.t, n_Cl = 1, d_dC = d.c, d_dE = d.e,
#                 b_Trt = FALSE, n_Seed = 1) # run for no treatment
# ResV_trt_Cpp =
#   MicroSimV_Cpp(v_S_t = v_M_1, t_P = t_p, v_C = c_vec, v_U = u_vec, n_I = n.i,
#                 n_S = n.s, n_T = n.t, n_Cl = 1, d_dC = d.c, d_dE = d.e,
#                 b_Trt = TRUE, n_Seed = 1) # run for treatment
#
# ###
#
# Default <- SS_MicroSim(n_i = 1000000)
SS_MicroSim(p_S1S2 = 0.155761, hr_S1 = 3.108213, hr_S2 = 6.030473)
###

library(devtools)
load_all()

samples_SS_MicroSim <- sample_prior_LHS(
  .n_samples = 100,
  .l_params = sickSicker_data$l_params)
samples1_SS_MicroSim <- sample_prior_FGS(
  .n_samples = 5,
  .l_params = sickSicker_data$l_params)
samples2_SS_MicroSim <- sample_prior_RGS(
  .n_samples = 50,
  .l_params = sickSicker_data$l_params)

GOF_wsse_SS_MicroSim <- wSSE_GOF(
  .func = SS_MicroSim, .optim = FALSE,
  .args = NULL,
  .samples = samples_SS_MicroSim,
  .l_targets = sickSicker_data$l_targets,
  .sample_method = "LHS")

GOF_llik_SS_MicroSim <- LLK_GOF(
  .func = SS_MicroSim, .optim = FALSE,
  .args = NULL,
  .samples = samples_SS_MicroSim,
  .l_targets = sickSicker_data$l_targets,
  .sample_method = "LHS")

samples <- sample_prior_LHS(
  .n_samples = 5,
  .l_params = sickSicker_data$l_params)

NM_optimise_wSSE_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000)

GB_optimise_wSSE_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000)

SA_optimise_wSSE_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

GA_optimise_wSSE_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'SSE',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

NM_optimise_LLK_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000)

GB_optimise_LLK_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000)

SA_optimise_LLK_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

GA_optimise_LLK_SS_MicroSim <- calibrateModel_directed(
  .l_params = sickSicker_data$l_params,
  .func = SS_MicroSim,
  .args = NULL,
  .gof = 'LLK',
  .samples = samples,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = sickSicker_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

samples <- sample_prior_LHS(
  .n_samples = 1000,
  .l_params = sickSicker_data$l_params)

SIR_SS_MicroSim = calibrateModel_beyesian(
  .b_method = 'SIR', .func = SS_MicroSim,
  .args = NULL,
  .l_targets = sickSicker_data$l_targets,
  .l_params = sickSicker_data$l_params, .samples = samples)

set.seed(1) # Function crashes on set.seed(1)
IMIS2_SS_MicroSim = calibrateModel_beyesian2(
  .b_method = 'IMIS', .func = SS_MicroSim,
  .args = NULL,
  .l_targets = sickSicker_data$l_targets,
  .l_params = sickSicker_data$l_params,
  .n_resample = 1000)

rm(likelihood, prior, sample.prior)
set.seed(1) # Function crashes on set.seed(1)
IMIS_SS_MicroSim = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = SS_MicroSim,
  .args = NULL,
  .l_targets = sickSicker_data$l_targets,
  .l_params = sickSicker_data$l_params,
  .n_resample = 1000)

l_optim_lists_SS_MicroSim <- list(
  GA_optimise_LLK_SS_MicroSim, GA_optimise_wSSE_SS_MicroSim,
  GB_optimise_LLK_SS_MicroSim, GB_optimise_wSSE_SS_MicroSim,
  NM_optimise_LLK_SS_MicroSim, NM_optimise_wSSE_SS_MicroSim,
  SA_optimise_LLK_SS_MicroSim, SA_optimise_wSSE_SS_MicroSim)

l_optim_lists_SS_MicroSim2 <- list(
  GOF_wsse_SS_MicroSim, GOF_llik_SS_MicroSim)

l_optim_lists_SS_MicroSim3 <- list(
  SIR_SS_MicroSim, IMIS_SS_MicroSim, IMIS2_SS_MicroSim)

PSA_values_SS_MicroSim <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_SS_MicroSim,
  .search_method = 'Directed',
  .PSA_runs = 1000,
  .l_params = sickSicker_data$l_params)

PSA_values_SS_MicroSim2 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_SS_MicroSim2,
  .search_method = 'Random',
  .PSA_runs = 1000,
  .l_params = sickSicker_data$l_params)

PSA_values_SS_MicroSim3 <- PSA_calib_values(
  .l_optim_lists = l_optim_lists_SS_MicroSim3,
  .search_method = 'Bayesian',
  .PSA_runs = 1000,
  .l_params = sickSicker_data$l_params)

#Models with transformations#############################################
## HID_markov
set.seed(1)
test <- HID_markov()
set.seed(1)
test2 <- HID_markov_2()
testthat::compare(test, test2)

betas <- rbeta(n = 1e+6, shape1 = 2, shape2 = 8)
prob_to_logit(mean(betas)); prob_to_logit(sd(betas))

logit_to_prob(prob_to_logit(mean(betas))); logit_to_prob(prob_to_logit(sd(betas)))
####
# Proper testing: #######################################################
library(devtools)
load_all()

# Transformed version of the HID_markov:
set.seed(1)
HID2_results <- list()
HID2_results$Prior_samples[['LHS']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data2$l_params)

HID2_results$Prior_samples[['FGS']] <- sample_prior_FGS(
  .n_samples = 10000,
  .l_params = HID_data2$l_params)

HID2_results$Prior_samples[['RGS']] <- sample_prior_RGS(
  .n_samples = 10000,
  .l_params = HID_data2$l_params)
##
HID2_results$Calib_results$Random[[1]] <- wSSE_GOF(
  .func = HID_markov_2, .optim = FALSE,
  .args = NULL,
  .samples = HID2_results$Prior_samples$LHS,
  .l_targets = HID_data2$l_targets,
  .sample_method = "LHS")

HID2_results$Calib_results$Random[[2]] <- LLK_GOF(
  .func = HID_markov_2, .optim = FALSE,
  .args = NULL,
  .samples = HID2_results$Prior_samples$LHS,
  .l_targets = HID_data2$l_targets,
  .sample_method = "LHS")
##
HID2_results$Prior_samples[['LHS_Directed']] <- sample_prior_LHS(
  .n_samples = 2,
  .l_params = HID_data2$l_params)

HID2_results$Calib_results$Directed[[1]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000)

HID2_results$Calib_results$Directed[[2]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000)

HID2_results$Calib_results$Directed[[3]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000)

HID2_results$Calib_results$Directed[[4]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000)

HID2_results$Calib_results$Directed[[5]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

HID2_results$Calib_results$Directed[[6]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

HID2_results$Calib_results$Directed[[7]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

HID2_results$Calib_results$Directed[[8]] <- calibrateModel_directed(
  .l_params = HID_data2$l_params,
  .func = HID_markov_2,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID2_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data2$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)
##
HID2_results$Prior_samples[['LHS_Bayesian']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data2$l_params)

HID2_results$Calib_results$Bayesian[[1]] = calibrateModel_beyesian(
  .b_method = 'SIR', .func = HID_markov_2,
  .args = NULL,
  .l_targets = HID_data2$l_targets,
  .n_resample = 10000,
  .l_params = HID_data2$l_params,
  .samples = HID2_results$Prior_samples$LHS_Bayesian)

set.seed(1) # Function crashes on set.seed(1)
HID2_results$Calib_results$Bayesian[[2]] = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = HID_markov_2,
  .args = NULL,
  .l_targets = HID_data2$l_targets,
  .l_params = HID_data2$l_params,
  .transform = TRUE,
  .n_resample = 10000,
  .IMIS_iterations = 400,
  .IMIS_sample = 100)
##
HID2_results$PSA_samples[["Directed"]] <- PSA_calib_values(
  .l_calib_res_lists = HID2_results$Calib_results$Directed,
  .search_method = 'Directed',
  .PSA_samples = 10000,
  .transform_ = TRUE,
  .l_params = HID_data2$l_params)

HID2_results$PSA_samples[["Random"]] <- PSA_calib_values(
  .l_calib_res_lists = HID2_results$Calib_results$Random,
  .search_method = 'Random',
  .PSA_samples = 10000,
  .transform_ = TRUE,
  .l_params = HID_data2$l_params)

HID2_results$PSA_samples[["Bayesian"]] <- PSA_calib_values(
  .l_calib_res_lists = HID2_results$Calib_results$Bayesian,
  .search_method = 'Bayesian',
  .PSA_samples = 10000,
  .transform_ = TRUE,
  .l_params = HID_data2$l_params)
##
HID2_results$PSA_results <- run_PSA(
  .func_ = HID_markov_2,
  .PSA_calib_values_ = c(HID2_results$PSA_samples$Directed,
                         HID2_results$PSA_samples$Random,
                         HID2_results$PSA_samples$Bayesian),
  .args_ = list(calibrate_ = FALSE,
                transform_ = FALSE),
  .PSA_unCalib_values_ = NULL)

# Untransformed version of HID_markov:
set.seed(1)
HID_results <- list()
HID_results$Prior_samples[['LHS']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)

HID_results$Prior_samples[['FGS']] <- sample_prior_FGS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)

HID_results$Prior_samples[['RGS']] <- sample_prior_RGS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)
##
HID_results$Calib_results$Random[[1]] <- wSSE_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = HID_results$Prior_samples$LHS[1:3,],
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")

HID_results$Calib_results$Random[[2]] <- LLK_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = HID_results$Prior_samples$LHS[1,],
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")
##
HID_results$Prior_samples[['LHS_Directed']] <- sample_prior_LHS(
  .n_samples = 5,
  .l_params = HID_data$l_params)

HID_results$Calib_results$Directed[[1]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

HID_results$Calib_results$Directed[[2]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'NM',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

HID_results$Calib_results$Directed[[3]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

HID_results$Calib_results$Directed[[4]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)

HID_results$Calib_results$Directed[[5]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)

HID_results$Calib_results$Directed[[6]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

HID_results$Calib_results$Directed[[7]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'LLK',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

HID_results$Calib_results$Directed[[8]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'SSE',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)
##
HID_results$Prior_samples[['LHS_Bayesian']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)

HID_results$Calib_results$Bayesian[[1]] = calibrateModel_beyesian(
  .b_method = 'SIR', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .n_resample = 10000,
  .l_params = HID_data$l_params,
  .samples = HID_results$Prior_samples$LHS_Bayesian)

set.seed(1) # Function crashes on set.seed(1)
HID_results$Calib_results$Bayesian[[2]] = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .l_params = HID_data$l_params,
  .transform = FALSE,
  .n_resample = 10000,
  .IMIS_iterations = 400,
  .IMIS_sample = 100)
##
HID_results$PSA_samples[["Directed"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Directed,
  .search_method = 'Directed',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)

HID_results$PSA_samples[["Random"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Random,
  .search_method = 'Random',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)

HID_results$PSA_samples[["Bayesian"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Bayesian,
  .search_method = 'Bayesian',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)
##
HID_results$PSA_results <- run_PSA(
  .func_ = HID_markov,
  .PSA_calib_values_ = c(HID_results$PSA_samples$Directed,
                         HID_results$PSA_samples$Random,
                         HID_results$PSA_samples$Bayesian),
  .args_ = list(calibrate_ = FALSE),
  .PSA_unCalib_values_ = NULL)

#Pete plot ################################
## install coda from CRAN and my MCIR package from github:
## devtools::install_github('petedodd/MCIR')
library(MCIR)
library(coda)

## rosebrock function log likelihood

## with gradient:
rosenbrock <- function(x){
  f <- (1-x[1])^2 + 100*(x[2] - x[1]^2)^2
  g <- c(100*2*(x[2] - x[1]^2)*(-2)*x[1] - 2*(1-x[1]),
         200*(x[2] - x[1]^2) )
  return(list(logp=-f,grad=-g))
}

## w/o gradient
rosen <- function(x) return(rosenbrock(x)$logp)

run <- nuts_da(rosengrad,5e3,1e3,runif(2))
corplot(run)

#R6: ###############################################################

cal = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2$l_params,
    .targets = HID_data2$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal$
  sampleR(
    .n_samples = 10000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal$
  calibrateR_random(
  .optim = FALSE,
  .maximise = TRUE,
  .weighted = TRUE,
  .sample_method = c("LHS", "RGS", "FGS"),
  .calibration_method = c("LLK", "SSE"))
cal$
  calibrateR_directed(
    .gof = c('LLK', 'SSE'),
    .n_samples = 5,
    .calibration_method = c('NM', 'BFGS', 'SANN', 'GA'),
    .sample_method = c("LHS", "RGS", "FGS"),
    .max_iterations = 1000,
    temp = 10,
    tmax = 10)
cal$
  calibrateR_bayesian(
  .b_method = c('SIR', 'IMIS'),
  .n_resample = 10000,
  .IMIS_iterations = 400,
  .IMIS_sample = 100)
cal$calibration_results$bayesian$SIR %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$SIR$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$calibration_results$bayesian$IMIS %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$IMIS$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$
  sample_PSA_values(
  .calibration_methods = c('Random', 'Directed', 'Bayesian'),
  .PSA_samples = 1000)
cal$
  run_PSA()
saveRDS(object = cal, file = "data/calibrateR_R6.rds")
cal <- readRDS(file = "data/calibrateR_R6.rds")
cal$
  calibrateR_directed(
    .gof = c('LLK'),
    .n_samples = 5,
    .calibration_method = c('BFGS'),
    .sample_method = c("RGS", "FGS"),
    .max_iterations = 1000,
    temp = 10,
    tmax = 10)

cal_HID_markov_2 = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2$l_params,
    .targets = HID_data2$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal_HID_markov_2$prior_samples = cal$prior_samples
cal_HID_markov_2$calibration_results = cal$calibration_results
cal_HID_markov_2$PSA_summary = cal$PSA_summary
cal_HID_markov_2$PSA_samples = cal$PSA_samples
cal_HID_markov_2$PSA_results = cal$PSA_results
cal_HID_markov_2$summarise_PSA()
cal_HID_markov_2$draw_plots()

# cal$
#   calibrateR_directed(
#     .gof = c('LLK'),
#     .n_samples = 1,
#     .calibration_method = c('GA', 'NM'),
#     .sample_method = c("LHS"),
#     .max_iterations = 1000,
#     temp = 10,
#     tmax = 10)

# testing:----
cal = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat$l_params,
    .targets = HID_data2_flat$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal$
  sampleR(
    .n_samples = 10000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
cal$calibration_results$bayesian$SIR %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$SIR$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$calibration_results$bayesian$IMIS %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$IMIS$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal$
  run_PSA()
cal$draw_plots()
saveRDS(object = cal, file = "data/calibrateR_R6_flat.rds")
# cal <- readRDS(file = "data/calibrateR_R6_flat.rds")
#
# cal_HID_markov_2 = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2$l_params,
#     .targets = HID_data2$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal_HID_markov_2$prior_samples = cal$prior_samples
# cal_HID_markov_2$calibration_results = cal$calibration_results
# cal_HID_markov_2$PSA_summary = cal$PSA_summary
# cal_HID_markov_2$PSA_samples = cal$PSA_samples
# cal_HID_markov_2$PSA_results = cal$PSA_results
# cal_HID_markov_2$summarise_PSA()
# cal$draw_plots()
# cal = cal_HID_markov_2
# exp(min(cal$prior_samples$LHS$mu_e))
# exp(max(cal$prior_samples$LHS$mu_e))
# exp(quantile(cal$prior_samples$LHS$mu_e, c(0.01, 0.99)))
# cal$prior_samples$LHS %>%
#   filter()
# exp(mean(cal$prior_samples$LHS$mu_e))
# exp(mean(cal$prior_samples$LHS$mu_e))

# testing 2 params:----
cal_2p = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat_2p$l_params,
    .targets = HID_data2_flat_2p$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal_2p$
  sampleR(
    .n_samples = 1000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal_2p$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 1000,
    temp = 10,
    tmax = 10)
cal_2p$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 1000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
cal_2p$calibration_results$bayesian$SIR %>%
  effective_sample_size(bayes_calib_output_list = .)
cal_2p$calibration_results$bayesian$SIR$Results %>%
  dplyr::distinct() %>%
  nrow()
cal_2p$calibration_results$bayesian$IMIS %>%
  effective_sample_size(bayes_calib_output_list = .)
cal_2p$calibration_results$bayesian$IMIS$Results %>%
  dplyr::distinct() %>%
  nrow()
cal_2p$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p$
  draw_plots()
cal_2p$
  run_PSA()
cal_2p$
  summarise_PSA()
#saveRDS(object = cal_2p, file = "data/calibrateR_R6_flat_2p.rds")
#cal_2p <- readRDS(file = "data/calibrateR_R6_flat_2p.rds")

# testing other 2 params:----
cal_2p2 = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat_2p2$l_params,
    .targets = HID_data2_flat_2p2$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal_2p2$
  sampleR(
    .n_samples = 1000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal_2p2$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p2$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 1000,
    temp = 10,
    tmax = 10)
cal_2p2$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 1000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
cal_2p2$calibration_results$bayesian$SIR %>%
  effective_sample_size(bayes_calib_output_list = .)
cal_2p2$calibration_results$bayesian$SIR$Results %>%
  dplyr::distinct() %>%
  nrow()
cal_2p2$calibration_results$bayesian$IMIS %>%
  effective_sample_size(bayes_calib_output_list = .)
cal_2p2$calibration_results$bayesian$IMIS$Results %>%
  dplyr::distinct() %>%
  nrow()
cal_2p2$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p2$draw_plots()
# cal_2p2$
#   run_PSA()
# saveRDS(object = cal_2p2, file = "data/calibrateR_R6_flat_2p2.rds")
# testing all ps:----
cal = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat$l_params,
    .targets = HID_data2_flat$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal$
  sampleR(
    .n_samples = 1000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 1000,
    temp = 10,
    tmax = 10)
cal$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 1000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
cal$calibration_results$bayesian$SIR %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$SIR$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$calibration_results$bayesian$IMIS %>%
  effective_sample_size(bayes_calib_output_list = .)
cal$calibration_results$bayesian$IMIS$Results %>%
  dplyr::distinct() %>%
  nrow()
cal$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal$
  run_PSA()
cal$draw_plots()
#saveRDS(object = cal, file = "data/calibrateR_R6_flat.rds")

## testing prior_posterior plots:----
# cal = readRDS(file = "inst/extdata/calibrateR_R6_flat_2p2.rds")
# cal_HID_markov_2 = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2_flat_2p2$l_params,
#     .targets = HID_data2_flat_2p2$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal_HID_markov_2$prior_samples = cal$prior_samples
# cal_HID_markov_2$calibration_results = cal$calibration_results
# cal_HID_markov_2$PSA_summary = cal$PSA_summary
# cal_HID_markov_2$PSA_samples = cal$PSA_samples
# cal_HID_markov_2$PSA_results = cal$PSA_results
# cal_HID_markov_2$draw_plots()
# cal_2p2 = cal_HID_markov_2
# saveRDS(object = cal_2p2, file = "inst/extdata/calibrateR_R6_flat_2p2.rds")

# cal = readRDS(file = "inst/extdata/calibrateR_R6_flat.rds")
# cal_HID_markov_2 = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2_flat$l_params,
#     .targets = HID_data2_flat$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal_HID_markov_2$prior_samples = cal$prior_samples
# cal_HID_markov_2$calibration_results = cal$calibration_results
# cal_HID_markov_2$PSA_summary = cal$PSA_summary
# cal_HID_markov_2$PSA_samples = cal$PSA_samples
# cal_HID_markov_2$PSA_results = cal$PSA_results
# cal_HID_markov_2$draw_plots()
# cal_HID_markov_2$print_test()
# cal = cal_HID_markov_2
# saveRDS(object = cal, file = "inst/extdata/calibrateR_R6_flat.rds")

# New testing 2 params:----
cal_2p = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat_2p$l_params,
    .targets = HID_data2_flat_2p$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal_2p$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal_2p$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 20,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal_2p$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal_2p$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal_2p$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal_2p$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p$
  draw_plots()
cal_2p$
  run_PSA()
cal_2p$
  summarise_PSA()
saveRDS(object = cal_2p, file = "inst/extdata/calibR_R6_flat_2p.rds")
#cal_2p <- readRDS(file = "data/calibrateR_R6_flat_2p.rds")

# New testing other 2 params:----
cal_2p2 = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat_2p2$l_params,
    .targets = HID_data2_flat_2p2$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal_2p2$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal_2p2$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p2$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 20,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal_2p2$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal_2p2$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p2$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal_2p2$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p2$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal_2p2$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p2$
  draw_plots()
cal_2p2$
  run_PSA()
cal_2p2$
  summarise_PSA()
saveRDS(object = cal_2p2, file = "inst/extdata/calibR_R6_flat_2p2.rds")
# New testing all ps:----
cal = calibR_R6$
  new(
    .model = HID_markov_2,
    .params = HID_data2_flat$l_params,
    .targets = HID_data2_flat$l_targets,
    .args = NULL,
    .transform = TRUE
  )
cal$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS", "RGS", "FGS")
  )
cal$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal$
  run_PSA()
cal$
  draw_plots()
cal$
  summarise_PSA()
saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat.rds")

# New testing 2 params (untransformed):----
cal_2p = calibR_R6$
  new(
    .model = HID_markov,
    .params = HID_data_flat_2p$l_params,
    .targets = HID_data_flat_2p$l_targets,
    .args = NULL,
    .transform = FALSE
  )
cal_2p$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS")
  )
cal_2p$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 20,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal_2p$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal_2p$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal_2p$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal_2p$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p$
  draw_plots()
cal_2p$
  run_PSA()
cal_2p$
  summarise_PSA()

saveRDS(object = cal_2p, file = "inst/extdata/calibR_R6_flat_2p_unT.rds")
rm(list = ls())
#cal_2p <- readRDS(file = "data/calibrateR_R6_flat_2p.rds")

# New testing other 2 params (untransformed):----
cal_2p2 = calibR_R6$
  new(
    .model = HID_markov,
    .params = HID_data_flat_2p2$l_params,
    .targets = HID_data_flat_2p2$l_targets,
    .args = NULL,
    .transform = FALSE
  )
cal_2p2$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS")
  )
cal_2p2$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal_2p2$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 20,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal_2p2$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal_2p2$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p2$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal_2p2$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal_2p2$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal_2p2$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal_2p2$
  draw_plots()
cal_2p2$
  run_PSA()
cal_2p2$
  summarise_PSA()

saveRDS(object = cal_2p2, file = "inst/extdata/calibR_R6_flat_2p2_unT.rds")
rm(list = ls())
# New testing all ps  (untransformed):----
cal = calibR_R6$
  new(
    .model = HID_markov,
    .params = HID_data_flat$l_params,
    .targets = HID_data_flat$l_targets,
    .args = NULL,
    .transform = FALSE
  )
cal$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS")
  )
cal$
  calibrateR_random(
    .optim = FALSE,
    .maximise = TRUE,
    .weighted = TRUE,
    .sample_method = "LHS",
    .calibration_method = "LLK")
cal$
  calibrateR_directed(
    .gof = 'LLK',
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10)
cal$
  calibrateR_bayesian(
    .b_method = c('SIR', 'IMIS'),
    .n_resample = 10000,
    .IMIS_iterations = 400,
    .IMIS_sample = 100)
# cal$calibration_results$bayesian$SIR %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal$calibration_results$bayesian$SIR$Results %>%
#   dplyr::distinct() %>%
#   nrow()
# cal$calibration_results$bayesian$IMIS %>%
#   effective_sample_size(bayes_calib_output_list = .)
# cal$calibration_results$bayesian$IMIS$Results %>%
#   dplyr::distinct() %>%
#   nrow()
cal$
  sample_PSA_values(
    .calibration_methods = c('Random', 'Directed', 'Bayesian'),
    .PSA_samples = 1000)
cal$
  run_PSA()
cal$
  draw_plots()
cal$
  summarise_PSA()

saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_unT.rds")
rm(list = ls())

##Amendments #################################################
# cal = calibR_R6$
#   new(
#     .model = HID_markov,
#     .params = HID_data_flat_2p$l_params,
#     .targets = HID_data_flat_2p$l_targets,
#     .args = NULL,
#     .transform = FALSE
#   )
# cal$prior_samples = calibR_R6_flat_2p_unT$prior_samples
# cal$calibration_results = calibR_R6_flat_2p_unT$calibration_results
# cal$PSA_samples = calibR_R6_flat_2p_unT$PSA_samples
# cal$PSA_results = calibR_R6_flat_2p_unT$PSA_results
# cal$PSA_summary = calibR_R6_flat_2p_unT$PSA_summary
# cal$draw_plots()
# cal$summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_2p_unT.rds")
# rm(list = ls())

# cal = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2_flat_2p$l_params,
#     .targets = HID_data2_flat_2p$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal$prior_samples = calibR_R6_flat_2p$prior_samples
# cal$calibration_results = calibR_R6_flat_2p$calibration_results
# cal$PSA_samples = calibR_R6_flat_2p$PSA_samples
# cal$PSA_results = calibR_R6_flat_2p$PSA_results
# cal$PSA_summary = calibR_R6_flat_2p$PSA_summary
# cal$draw_plots()
# cal$summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_2p.rds")
# rm(list = ls())

# cal = calibR_R6$
#   new(
#     .model = HID_markov,
#     .params = HID_data_flat_2p2$l_params,
#     .targets = HID_data_flat_2p2$l_targets,
#     .args = NULL,
#     .transform = FALSE
#   )
# cal$prior_samples = calibR_R6_flat_2p2_unT$prior_samples
# cal$calibration_results = calibR_R6_flat_2p2_unT$calibration_results
# cal$PSA_samples = calibR_R6_flat_2p2_unT$PSA_samples
# cal$PSA_results = calibR_R6_flat_2p2_unT$PSA_results
# cal$PSA_summary = calibR_R6_flat_2p2_unT$PSA_summary
# cal$draw_plots()
# cal$summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_2p2_unT.rds")
# rm(list = ls())

# cal = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2_flat_2p2$l_params,
#     .targets = HID_data2_flat_2p2$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal$prior_samples = calibR_R6_flat_2p2$prior_samples
# cal$calibration_results = calibR_R6_flat_2p2$calibration_results
# cal$PSA_samples = calibR_R6_flat_2p2$PSA_samples
# cal$PSA_results = calibR_R6_flat_2p2$PSA_results
# cal$PSA_summary = calibR_R6_flat_2p2$PSA_summary
# cal$draw_plots()
# cal$summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_2p2.rds")
# rm(list = ls())

# cal = calibR_R6$
#   new(
#     .model = HID_markov,
#     .params = HID_data_flat$l_params,
#     .targets = HID_data_flat$l_targets,
#     .args = NULL,
#     .transform = FALSE
#   )
# cal$prior_samples = calibR_R6_flat_unT$prior_samples
# cal$calibration_results = calibR_R6_flat_unT$calibration_results
# cal$PSA_samples = calibR_R6_flat_unT$PSA_samples
# cal$PSA_results = calibR_R6_flat_unT$PSA_results
# cal$PSA_summary = calibR_R6_flat_unT$PSA_summary
# cal$
#   calibrateR_bayesian(
#     .b_method = c('SIR', 'IMIS'),
#     .n_resample = 10000,
#     .IMIS_iterations = 400,
#     .IMIS_sample = 100)
# cal$
#   sample_PSA_values(
#     .calibration_methods = c('Random', 'Directed', 'Bayesian'),
#     .PSA_samples = 1000)
# cal$
#   run_PSA()
# cal$
#   draw_plots()
# cal$
#   summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat_unT.rds")
# rm(list = ls())

# cal = calibR_R6$
#   new(
#     .model = HID_markov_2,
#     .params = HID_data2_flat$l_params,
#     .targets = HID_data2_flat$l_targets,
#     .args = NULL,
#     .transform = TRUE
#   )
# cal$prior_samples = calibR_R6_flat$prior_samples
# cal$calibration_results = calibR_R6_flat$calibration_results
# cal$PSA_samples = calibR_R6_flat$PSA_samples
# cal$PSA_results = calibR_R6_flat$PSA_results
# cal$PSA_summary = calibR_R6_flat$PSA_summary
# cal$draw_plots()
# cal$summarise_PSA()
#
# saveRDS(object = cal, file = "inst/extdata/calibR_R6_flat.rds")
# rm(list = ls())


#Rosenbrok################################

cal <- calibR_R6$new(
  .model = calibTest_rosen,
  .args = NULL,
  .params = ROSEN_data$l_params,
  .targets = ROSEN_data$l_targets,
  .transform = FALSE)

cal$
  sampleR(
    .n_samples = 100000,
    .sampling_method = c("LHS")
  )
cal$
  calibrateR_directed(
    .gof = 'Rosen',
    .gof_func = calibTest_rosen,
    .n_samples = 10,
    .calibration_method = c('NM', 'BFGS', 'SANN'),
    .sample_method = "LHS",
    .max_iterations = 10000,
    temp = 10,
    tmax = 10,
    .maximise = FALSE)