R/run_sim.R

In WARDEN: Workflows for Health Technology Assessments in R using Discrete EveNts

#' Run the simulation
#'
#' @param arm_list A vector of the names of the interventions evaluated in the simulation
#' @param sensitivity_inputs A list of sensitivity inputs that do not change within a sensitivity in a similar fashion to common_all_inputs, etc
#' @param common_all_inputs A list of inputs common across patients that do not change within a simulation
#' @param common_pt_inputs A list of inputs that change across patients but are not affected by the intervention
#' @param unique_pt_inputs A list of inputs that change across each intervention
#' @param init_event_list A list of initial events and event times. If no initial events are given, a "Start" event at time 0 is created automatically
#' @param evt_react_list A list of event reactions
#' @param util_ongoing_list Vector of QALY named variables that are accrued at an ongoing basis (discounted using drq)
#' @param util_instant_list Vector of QALY named variables that are accrued instantaneously at an event (discounted using drq)
#' @param util_cycle_list Vector of QALY named variables that are accrued in cycles (discounted using drq)
#' @param cost_ongoing_list Vector of cost named variables that are accrued at an ongoing basis (discounted using drc)
#' @param cost_instant_list Vector of cost named variables that are accrued instantaneously at an event (discounted using drc)
#' @param cost_cycle_list Vector of cost named variables that are accrued in cycles (discounted using drc)
#' @param other_ongoing_list Vector of other named variables that are accrued at an ongoing basis (discounted using drq)
#' @param other_instant_list Vector of other named variables that are accrued instantaneously at an event (discounted using drq)
#' @param npats The number of patients to be simulated (it will simulate npats * length(arm_list))
#' @param n_sim The number of simulations to run per sensitivity
#' @param psa_bool A boolean to determine if PSA should be conducted. If n_sim > 1 and psa_bool = FALSE, the differences between simulations will be due to sampling
#' @param sensitivity_bool A boolean to determine if Scenarios/DSA should be conducted. 
#' @param sensitivity_names A vector of scenario/DSA names that can be used to select the right sensitivity (e.g., c("Scenario_1", "Scenario_2")). The parameter "sens_name_used" is created from it which corresponds to the one being used for each iteration.
#' @param n_sensitivity Number of sensitivity analysis (DSA or Scenarios) to run. It will be interacted with sensitivity_names argument if not null (n_sensitivityitivity = n_sensitivity * length(sensitivity_names)). For DSA, it should be as many parameters as there are. For scenario, it should be 1.
#' @param input_out A vector of variables to be returned in the output data frame
#' @param ipd Integer taking value 1 for full IPD data returned, and 2 IPD data but aggregating events (returning last value for numeric/character/factor variables. For other objects (e.g., matrices), the IPD will still be returned as the aggregation rule is not clear). Other values mean no IPD data returned (removes non-numerical or length>1 items)
#' @param timed_freq If NULL, it does not produce any timed outputs. Otherwise should be a number (e.g., every 1 year)
#' @param debug If TRUE, will generate a log file
#' @param accum_backwards If TRUE, the ongoing accumulators will count backwards (i.e., the current value is applied until the previous update). If FALSE, the current value is applied between the current event and the next time it is updated.
#' @param continue_on_error If TRUE, on error it will attempt to continue by skipping the current simulation 
#' @param seed Starting seed to be used for the whole analysis. If null, it's set to 1 by default.
#'
#' @return A list of data frames with the simulation results
#' @importFrom progressr with_progress
#' @importFrom progressr handlers
#' @importFrom progressr handler_txtprogressbar
#' @importFrom progressr progressor
#'
#' @export
#' @details This function is slightly different from `run_sim_parallel`.
#' `run_sim_parallel` only runs multiple-core at the simulation level.
#' `run_sim` uses only-single core.
#' `run_sim` can be more efficient if using only one simulation (e.g., deterministic),
#'  while `run_sim_parallel` will be more efficient if the number of simulations is >1 (e.g., PSA).
#'  
#'  Event ties are processed in the order declared within the `init_event_list` argument (`evts` argument within the first sublist of that object).
#'  To do so, the program automatically adds a sequence from to 0 to the (number of events - 1) times 1e-10 to add to the event times when selecting the event with minimum time.
#'  This time has been selected as it's relatively small yet not so small as to be ignored by which.min (see .Machine for more details)
#'  
#'  A list of protected objects that should not be used by the user as input names  or in the global environment to avoid the risk of overwriting them is as follows:
#'  c("arm", "arm_list", "categories_for_export", "cur_evtlist", "curtime", "evt", "i", "prevtime", "sens", "simulation", "sens_name_used","list_env","uc_lists","npats","ipd").
#'  
#'  The engine uses the L'Ecuyer-CMRG for the random number generator. 
#'  Note that the random seeds are set to be unique in their category (i.e., at patient level, patient-arm level, etc.)
#'  
#'  If no `drc` or `drq parameters are passed within any of the input lists, these are assigned value 0.03.
#'  
#'  Ongoing items will look backward to the last time updated when performing the discounting and accumulation.
#'  This means that the user does not necessarily need to keep updating the value, but only add it when the value
#'  changes looking forward (e.g., o_q = utility at event 1, at event 2 utility does not change, but at event 3 it does,
#'  so we want to make sure to add o_q = utility at event 3 before updating utility. The program will automatically 
#'  look back until event 1). Note that in previous versions of the package backward was the default, and now this has switched to forward.
#'  
#'  It is important to note that the QALYs and Costs (ongoing or instant or per cycle) used should be of length 1. 
#'  If they were of length > 1, the model would expand the data,
#'   so instead of having each event as a row, the event would have N rows (equal to the length of the costs/qalys to discount passed). 
#'   This means more processing of the results data would be needed in order for it to provide the correct results.
#'  
#'  If the `cycle` lists are used, then it is expected the user will declare as well the name of the variable
#'   pasted with `cycle_l` and `cycle_starttime` (e.g., c_default_cycle_l and c_default_cycle_starttime) to 
#'   ensure the discounting can be computed using cycles, with cycle_l being the cycle length, and cycle_starttime 
#'   being the starting time in which the variable started counting.
#'   
#'  `debug = TRUE` will export a log file with the timestamp up the error in the main working directory.
#'  
#'  `continue_on_error` will skip the current simulation (so it won't continue for the rest of patient-arms) if TRUE.
#'   Note that this will make the progress bar not correct, as a set of patients that were expected to be run is not.
#'
#' @examples
#' library(magrittr)
#' common_all_inputs <-add_item(
#' util.sick = 0.8,
#' util.sicker = 0.5,
#' cost.sick = 3000,
#' cost.sicker = 7000,
#' cost.int = 1000,
#' coef_noint = log(0.2),
#' HR_int = 0.8,
#' drc = 0.035, #different values than what's assumed by default
#' drq = 0.035,
#' random_seed_sicker_i = sample.int(100000,5,replace = FALSE)
#' )
#' 
#' common_pt_inputs <- add_item(death= max(0.0000001,rnorm(n=1, mean=12, sd=3))) 
#' 
#' unique_pt_inputs <- add_item(fl.sick = 1,
#'                              q_default = util.sick,
#'                              c_default = cost.sick + if(arm=="int"){cost.int}else{0}) 
#'                              
#' init_event_list <- 
#' add_tte(arm=c("noint","int"), evts = c("sick","sicker","death") ,input={
#'   sick <- 0
#'   sicker <- draw_tte(1,dist="exp",
#'    coef1=coef_noint, beta_tx = ifelse(arm=="int",HR_int,1),
#'     seed = random_seed_sicker_i[i])
#'   
#' })   
#' 
#' evt_react_list <-
#' add_reactevt(name_evt = "sick",
#'              input = {}) %>%
#'   add_reactevt(name_evt = "sicker",
#'                input = {
#'                  modify_item(list(q_default = util.sicker,
#'                                   c_default = cost.sicker + if(arm=="int"){cost.int}else{0},
#'                                   fl.sick = 0)) 
#'                }) %>%
#'   add_reactevt(name_evt = "death",
#'                input = {
#'                  modify_item(list(q_default = 0,
#'                                   c_default = 0, 
#'                                   curtime = Inf)) 
#'                }) 
#'                
#' util_ongoing <- "q_default"
#' cost_ongoing <- "c_default"
#'                           
#' 
#' run_sim(arm_list=c("int","noint"),
#' common_all_inputs = common_all_inputs,
#' common_pt_inputs = common_pt_inputs,
#' unique_pt_inputs = unique_pt_inputs,
#' init_event_list = init_event_list,
#' evt_react_list = evt_react_list,
#' util_ongoing_list = util_ongoing,
#' cost_ongoing_list = cost_ongoing,
#' npats = 2,
#' n_sim = 1,
#' psa_bool = FALSE,
#' ipd = 1)
#' 

run_sim <- function(arm_list=c("int","noint"),
                   sensitivity_inputs=NULL,
                   common_all_inputs=NULL,
                   common_pt_inputs=NULL,
                   unique_pt_inputs=NULL,
                   init_event_list = NULL,
                   evt_react_list = evt_react_list,
                   util_ongoing_list = NULL,
                   util_instant_list = NULL,
                   util_cycle_list = NULL,
                   cost_ongoing_list = NULL,
                   cost_instant_list = NULL,
                   cost_cycle_list = NULL,
                   other_ongoing_list = NULL,
                   other_instant_list = NULL,
                   npats=500,
                   n_sim=1,
                   psa_bool = NULL,
                   sensitivity_bool = FALSE,
                   sensitivity_names = NULL,
                   n_sensitivity = 1,
                   input_out = NULL,
                   ipd = 1,
                   timed_freq = NULL,
                   debug = FALSE,
                   accum_backwards = FALSE,
                   continue_on_error = FALSE,
                   seed = NULL){



# Set-up basics -----------------------------------------------------------
  
  #Store original rng kind and use L'Ecuyer-CMRG 
  rng_kind_store <- RNGkind()[1]
  RNGkind("L'Ecuyer-CMRG")
  
  #Stop simulation if forbidden objects are used.
  list_forbidden_names <-  c("arm", "arm_list", "categories_for_export", "cur_evtlist", "curtime", "evt", "i", "prevtime", "sens", "simulation", "sens_name_used","list_env","uc_lists","npats","ipd") 
  matched_list_forbidden <- which(list_forbidden_names %in% c(
    ls(parent.frame()),
    names(common_all_inputs),
    names(sensitivity_inputs),
    names(common_pt_inputs),
    names(unique_pt_inputs),
    names(arm_list),
    names(evt_react_list)
  )
  )
  
  if(  length(matched_list_forbidden)>0){
    stop(paste0("Name(s) or object: `", list_forbidden_names[matched_list_forbidden],"` defined belong to the list of forbidden names, which can cause issues in the model.
    This is likely caused by an object being declared in the parent/global environment (e.g., `i`) or in the inputs. 
    You can use functions like find('i') to help you locate them. 
    Please remove or rename those objects. See run_sim or run_sim_parallel for a full list of these names.\n"))
  }
  
  if(length(names(evt_react_list)) != length(init_event_list[[1]][["evts"]])){
    dif <- setdiff(names(evt_react_list),init_event_list[[1]][["evts"]])
    stop(paste0("Number of events defined in init_event_list first sublist `evts` is not equal to the number of events defined in evt_react_list.
         Please ensure both have equal length, ",dif," is/are missing from one of the lists."))
  }
  
  arm_list <- arm_list #this is done as otherwise there are problems passing arguments from function to function
  
  #get cost/utility categories
  categories_costs_ongoing <- cost_ongoing_list
  categories_costs_instant <- cost_instant_list
  categories_costs_cycle   <- cost_cycle_list

  categories_utilities_ongoing <- util_ongoing_list
  categories_utilities_instant <- util_instant_list
  categories_utilities_cycle   <- util_cycle_list
 
  categories_other_ongoing <- other_ongoing_list
  categories_other_instant <- other_instant_list
  
  
  #Remove NULL values
  categories_for_export <- unique(
    c(if(!is.null(categories_costs_ongoing)){categories_costs_ongoing},
      if(!is.null(categories_costs_instant)){categories_costs_instant},
      if(!is.null(categories_costs_cycle)){categories_costs_cycle},
      if(!is.null(categories_utilities_ongoing)){categories_utilities_ongoing},
      if(!is.null(categories_utilities_instant)){categories_utilities_instant},
      if(!is.null(categories_utilities_cycle)){categories_utilities_cycle},
      if(!is.null(categories_other_ongoing)){categories_other_ongoing},
      if(!is.null(categories_other_instant)){categories_other_instant}
    ))
  
  output_sim <- list()
  
  final_log <- list()
  log_list <- list()
  .skip_to_next <- FALSE
  flag_noerr_sim <- TRUE
  
  start_time <-  proc.time()
  
# Analysis loop ---------------------------------------------------------

  if (is.null(sensitivity_names)) {
    length_sensitivities <- n_sensitivity
  } else{
    length_sensitivities <- n_sensitivity * length(sensitivity_names)
  }
    
    progressr::handlers(progressr::handler_txtprogressbar(width=100))
    
    #Progress is forced here, because otherwise the user could not see this evolve properly.
    progressr::with_progress({
     pb <- progressr::progressor(min(npats*length_sensitivities*n_sim,50)) 
      
  #Need to figure out how to distinguish DSA (as many sensitivities as parameters) and Scenarios (as many sensivities as scenarios)
  for (sens in 1:length_sensitivities) {
    message(paste0("Analysis number: ",sens))
    
    tryCatch({
    
    start_time_analysis <-  proc.time()

    output_sim[[sens]] <- list() #initialize sensitivity lists
    
    #e.g., if length_sensitivities is 50 (25 param x 2 DSAs) then take at each iteration the divisor to see which name should be applied
    if (!is.null(sensitivity_names)) {
    sens_name_used <- sensitivity_names[ceiling(sens/n_sensitivity)] 
    } else{
    sens_name_used <- ""
    }
    
    input_list_sens <- list(
                       psa_bool = psa_bool,
                       init_event_list = init_event_list,
                       precision_times = if(!is.null(init_event_list)){setNames((seq_along(init_event_list[[1]][["evts"]])-1)*1e-10,
                                                                                init_event_list[[1]][["evts"]])},
                       evt_react_list = evt_react_list,
                       uc_lists = list(util_ongoing_list = util_ongoing_list,
                                       util_instant_list = util_instant_list,
                                       util_cycle_list = util_cycle_list,
                                       cost_ongoing_list = cost_ongoing_list,
                                       cost_instant_list = cost_instant_list,
                                       cost_cycle_list = cost_cycle_list,
                                       other_ongoing_list = other_ongoing_list,
                                       other_instant_list = other_instant_list,
                                       ongoing_inputs = c(
                                         if(!is.null(categories_costs_ongoing)){categories_costs_ongoing},
                                         if(!is.null(categories_utilities_ongoing)){categories_utilities_ongoing},
                                         if(!is.null(categories_other_ongoing)){categories_other_ongoing}
                                       ),
                                       instant_inputs = c(
                                         if(!is.null(categories_costs_instant)){categories_costs_instant},
                                         if(!is.null(categories_utilities_instant)){categories_utilities_instant},
                                         if(!is.null(categories_other_instant)){categories_other_instant}
                                       ),
                                       cycle_inputs = c(
                                         if(!is.null(categories_costs_cycle)){categories_costs_cycle},
                                         if(!is.null(categories_utilities_cycle)){categories_utilities_cycle}
                                       ),
                                       cost_categories_ongoing = categories_costs_ongoing,
                                       l_cost_categories_ongoing = length(categories_costs_ongoing),
                                       cost_categories_instant = categories_costs_instant,
                                       l_cost_categories_instant = length(categories_costs_instant),
                                       cost_categories_cycle = categories_costs_cycle,
                                       l_cost_categories_cycle = length(categories_costs_cycle),
                                       util_categories_ongoing = categories_utilities_ongoing,
                                       l_util_categories_ongoing = length(categories_utilities_ongoing),
                                       util_categories_instant = categories_utilities_instant,
                                       l_util_categories_instant = length(categories_utilities_instant),
                                       util_categories_cycle = categories_utilities_cycle,
                                       l_util_categories_cycle = length(categories_utilities_cycle),
                                       other_categories_ongoing = categories_other_ongoing,
                                       l_other_categories_ongoing = length(categories_other_ongoing),
                                       other_categories_instant = categories_other_instant,
                                       l_other_categories_instant = length(categories_other_instant)
                       ),
                       input_out = unique(c(input_out,categories_for_export)),
                       categories_for_export = categories_for_export,
                       ipd = ipd,
                       arm_list = arm_list,
                       npats = npats,
                       n_sim = n_sim,
                       sensitivity_bool = sensitivity_bool,
                       n_sensitivity = n_sensitivity,
                       length_sensitivities = length_sensitivities,
                       sens = sens,
                       sensitivity_names = sensitivity_names,
                       sens_name_used = sens_name_used,
                       timed_freq = timed_freq,
                       debug = debug,
                       accum_backwards = accum_backwards,
                       continue_on_error = continue_on_error,
                       log_list = list()
                      )
    
    if(is.null(seed)){
      seed <- 1
    }
    set.seed(seed)
    
    
    # Draw Common parameters  -------------------------------
    if(!is.null(sensitivity_inputs)){
      
      input_list_sens <- load_inputs(inputs = input_list_sens,list_uneval_inputs = sensitivity_inputs)
      
      
      if(input_list_sens$debug){ 
        names_sens_input <- names(sensitivity_inputs)
        prev_value <- setNames(vector("list", length(sensitivity_inputs)), names_sens_input)
        dump_info <- list(
          list(
            prev_value = prev_value,
            cur_value  = input_list_sens[names_sens_input]
          )
        )
        
        names(dump_info) <- paste0("Analysis: ", input_list_sens$sens," ", input_list_sens$sens_name_used,
                                   "; Structural"
        )
        
        log_list <- c(log_list,dump_info)
      }
      
    }
    
    #Make sure there are no duplicated inputs in the model, if so, take the last one
    duplic <- duplicated(names(input_list_sens),fromLast = TRUE)
    if (sum(duplic)>0 & sens==1) { warning("Duplicated items detected in the Analysis, using the last one added.\n")  }
    input_list_sens <- input_list_sens[!duplic]

# Simulation loop ---------------------------------------------------------

      
    for (simulation in 1:n_sim) {
      message(paste0("Simulation number: ",simulation))
      
      tryCatch({
        
      
      start_time_sim <-  proc.time()
      
      input_list <- c(input_list_sens,list(simulation=simulation))
      
      set.seed(simulation*1007*seed)
      
      # Draw Common parameters  -------------------------------
      if(!is.null(common_all_inputs)){
        
        input_list <- load_inputs(inputs = input_list,list_uneval_inputs = common_all_inputs)

        if(input_list_sens$debug){ 
          names_all_input <- names(common_all_inputs)
          prev_value <- setNames(vector("list", length(common_all_inputs)), names_all_input)
          prev_value[names_all_input] <- input_list_sens[names_all_input]
          dump_info <- list(
            list(
              prev_value = prev_value,
              cur_value  = input_list[names_all_input]
            )
          )
          
          names(dump_info) <- paste0("Analysis: ", input_list$sens," ", input_list$sens_name_used,
                                     "; Sim: ", input_list$sim,
                                     "; Statics"
          )
          
          log_list <- c(log_list,dump_info)
        }
      }
      
      #Make sure there are no duplicated inputs in the model, if so, take the last one
      duplic <- duplicated(names(input_list),fromLast = TRUE)
      if (sum(duplic)>0 & simulation==1 & sens==1) { warning("Duplicated items detected in the Simulation, using the last one added.\n")  }
      input_list <- input_list[!duplic]
      
      if(is.null(input_list$drc)){input_list$drc <- 0.03}
      if(is.null(input_list$drq)){input_list$drq <- 0.03}
      
      # Run engine ----------------------------------------------------------
  
        final_output <- run_engine(arm_list=arm_list,
                                        common_pt_inputs=common_pt_inputs,
                                        unique_pt_inputs=unique_pt_inputs,
                                        input_list = input_list,
                                        pb = pb,
                                        seed = seed)                    
      
      if(!is.null(final_output$error_m)){
        if((n_sim > 1 | n_sensitivity > 1) & continue_on_error){
          message(if(debug){"Log will be exported. "},
                  "Continuing on error, error message at analysis ",
                  sens,
                  "; simulation: ",
                  if(exists("simulation")){simulation}else{"None"},
                  ". Error message: ",final_output$error_m
                  )
          next
        } else{
        stop(final_output$error_m)
        }
      } 
      if (input_list$ipd>0) {
        final_output$merged_df$simulation <- simulation
        final_output$merged_df$sensitivity <- sens
      }
      
      final_output <- c(list(sensitivity_name = sens_name_used), final_output)
      
      
      
      if(debug){
        log_list <- lapply(log_list,transform_debug)
        
        final_output$log_list <- c(log_list,final_output$log_list)
      }
      
      output_sim[[sens]][[simulation]] <- final_output
  
      message(paste0("Time to run simulation ", simulation,": ",  round(proc.time()[3]- start_time_sim[3] , 2 ), "s"))
      
      
      }, error = function(e) {
        if(continue_on_error){
          .skip_to_next <<- TRUE
          message(if(debug){"Log will be exported. "},
                  "Continuing on error, error message at analysis ",
                  sens,
                  "; simulation: ",
                  if(exists("simulation")){simulation}else{"None"},
                  ". Error message: ",e$message)
        } else{
          if(debug){
            flag_noerr_sim <<- FALSE
            
            if(!exists("final_output")){
              export_log(lapply(log_list,transform_debug),paste0("log_model_",format(Sys.time(), "%Y_%m_%d_%Hh_%mm_%Ss"),".txt"))
              stop(e$message)
            }else{
              if(length(log_list)>0){
                log_list <- lapply(log_list,transform_debug)
                final_output <- list()
                final_output$log_list <- c(log_list,final_output$log_list)
                
                if(exists("output_sim")){
                  if(!exists("simulation")){
                    simulation <- 1
                  }
                  output_sim[[sens]][[simulation]] <- final_output
                  final_log <- unlist(
                    unlist(
                      lapply(output_sim, function(y) lapply(y, function(x) x$log_list )),
                      recursive = FALSE),
                    recursive = FALSE)
                } else{
                  final_log <- final_output$log_list
                }
                
                export_log(final_log,paste0("log_model_",format(Sys.time(), "%Y_%m_%d_%Hh_%mm_%Ss"),".txt"))
              }
            }
            stop("Log Exported. Error message at analysis ",
                  sens,
                  "; simulation: ",
                  if(exists("simulation")){simulation}else{"None"},
                  ". Error message: ",
                  e$message)
          }else{
            stop("Error message at analysis ",
                 sens,
                 "; simulation: ",
                 if(exists("simulation")){simulation}else{"None"},
                 ". Error message: ",
                 e$message)
          }
        }
      }
      )
      if(.skip_to_next) { next } 
      
    }

    message(paste0("Time to run analysis ", sens,": ",  round(proc.time()[3]- start_time_analysis[3] , 2 ), "s"))
    
    
    }, error = function(e) {
      if(continue_on_error){
        .skip_to_next <<- TRUE
        message(if(debug){"Log will be exported. "},
                "Continuing on error, error message at analysis ",
                sens,
                "; simulation: ",
                if(exists("simulation")){simulation}else{"None"},
                ". Error message: ",
                e$message)
        } else{
        if(debug & flag_noerr_sim){
          if(length(log_list)>0){
            if(!exists("final_output")){
              export_log(lapply(log_list,transform_debug),paste0("log_model_",format(Sys.time(), "%Y_%m_%d_%Hh_%mm_%Ss"),".txt"))
              stop(e$message)
            }else{
                  if(exists("output_sim")){
                    final_log <- unlist(
                      unlist(
                        lapply(output_sim, function(y) lapply(y, function(x) x$log_list )),
                        recursive = FALSE),
                      recursive = FALSE)
                  } else{
                    log_list <- lapply(log_list,transform_debug)
                    final_output <- list()
                    final_output$log_list <- c(log_list,final_output$log_list)
                    final_log <- final_output$log_list
                  }
                  
                  export_log(final_log,paste0("log_model_",format(Sys.time(), "%Y_%m_%d_%Hh_%mm_%Ss"),".txt"))
              } 
          } else{
            message("No data to export.")
          }
          stop("Log will be exported. Error message at analysis ",
               sens,
               "; simulation: ",
               if(exists("simulation")){simulation}else{"None"},
               ". Error message: ",
               e$message)
        }else{
          stop(e$message)
        }
      }
    }
  )
    if(.skip_to_next) { next } 
    
    }
    
     message(paste0("Total time to run: ",  round(proc.time()[3]- start_time[3] , 2), "s"))


  # Export results ----------------------------------------------------------
  if(debug){
    if(length(log_list)>0){
      if(exists("output_sim")){
        final_log <- unlist(
          unlist(
            lapply(output_sim, function(y) lapply(y, function(x) x$log_list )),
            recursive = FALSE),
          recursive = FALSE)
      } else{
        log_list <- lapply(log_list,transform_debug)
        final_output <- list()
        final_output$log_list <- c(log_list,final_output$log_list)
        final_log <- final_output$log_list
      }
    
    export_log(final_log,paste0("log_model_",format(Sys.time(), "%Y_%m_%d_%Hh_%mm_%Ss"),".txt"))
    } else{
      message("No data to export.")
    }
  }
  
  
  results <- output_sim
  
  #Retore original rng kind
  RNGkind(rng_kind_store)
  
   }, enable=TRUE)
    
  
  return(results)
  

}