R/simulator.R

In yuehmeir2/CATSimulator: Library for Simulation of Adaptive Item Selection

#' Return a generator function that will generate one simulee each time it is called.
#' The true theta of the generated simulees will follow a normal distribution.
#' Once the specified number of simulees has been reached, no more simulees will be generated.
#'
#' @param numSimulees The number of simulees that should be generated.
#' @param mean The mean of the true theta for the generated simulees.
#' @param sd The standard deviation of the true theta for the generated simulees.
#' @param randomSeeds A vector of integers. The source of random seeds for each simulee.
#' @return A generator function that will generate one simulee each time it is called.
#' @examples
#'   simulees = generateSimuleesByNormal(numSimulees = 5, mean = 0.5, sd = 1.2, randomSeeds = 10001:10003)
#' @export
generateSimuleesByNormal <- function (numSimulees, mean = 0, sd = 1, randomSeeds = NULL, min_value = -Inf, max_value = Inf) {
  if (!is.null(randomSeeds)) set.seed(randomSeeds[[1]])
  trueThetas = rnorm(numSimulees, mean, sd) %>% pmin(max_value) %>% pmax(min_value)

  simulee = generateSimuleesByTrueTheta(trueThetas, randomSeeds)

  return(simulee)
}

#' Return a generator function that will generate one simulee each time it is called.
#' The true theta of the generated simulees will come from the next value in the given vector.
#' Once all true thetas have been consumed, no more simulees will be generated.
#'
#' @param trueThetas A vector of doubles. The source of true thetas for each simulee.
#' @param randomSeeds A vector of integers. The source of random seeds for each simulee.
#' @return A generator function that will generate one simulee each time it is called.
#' @examples
#'   simulees = generateSimuleesByTrueTheta(-2:2, 10001:10003)
#' @import tidyverse
#' @export
generateSimuleesByTrueTheta <- function (trueThetas, randomSeeds = NULL) {
  suppressPackageStartupMessages(require(tidyverse))

  # generate simulee ids with leading zeroes required to fit the number of trueThetas
  id_pattern = paste0("sim%0", ceiling(log10(length(trueThetas)+1)), "d")
  id = sprintf(id_pattern, 1:length(trueThetas))

  if (is.null(randomSeeds) || (length(randomSeeds) < 1)) {
    # generate randomSeeds randomly
    randomSeeds = runif(length(trueThetas), 10000, .Machine$integer.max - 3000)
  } else if (length(randomSeeds) != length(trueThetas)) {
    # recycle the randomSeeds values to match the length of trueThetas
    randomSeeds = rep(randomSeeds, ceiling(length(trueThetas) / length(randomSeeds)))[1:length(trueThetas)]
  }

  simulee = tibble(
    ID = as.character(id),
    SEED = as.integer(randomSeeds),
    TRUE_THETA = as.double(trueThetas)
  )

  return(simulee)
}

#' initSimulation
#'
#' @examples
#'   options(warn=-1)
#'   simulation = readRDS(system.file("example/passage-optimal.rds", package = "CATSimulator"))
#'   simulation = initSimulation(simulation)
#' @import tidyverse
#' @export
initSimulation <- function (simulation) {
  suppressPackageStartupMessages(require(tidyverse))
  # simulation = readRDS("/Users/ychien/Documents/code/R/MSTPrototype/MSTSuite/inst/example/mst12_p1_cons6.rds")
  # Generate 'item_pars' matrix including only the IRT params portion of the itempool
  if(!"item_pars" %in% names(simulation)) {
    simulation$item_pars = as.matrix(simulation$itempool[,startsWith(colnames(simulation$itempool), "PAR_")])
    dimnames(simulation$item_pars) = list(simulation$itempool$ITEM_ID, colnames(simulation$itempool)[startsWith(colnames(simulation$itempool), "PAR_")])
  }

  # Generate modules
  if ((!"modules" %in% names(simulation)) || (("generateModules" %in% names(simulation$control)) && simulation$control$generateModules)) {
    # Create an empty "modules" tibble, if necessary
    if (!"modules" %in% names(simulation)) {
      simulation$modules = tibble(
        MODULE_ID = character(),
        ITEM_IDS = list()
      )
    }

    # Generate a row for each item which does not already belong to a module.
    orphanItem = rep(TRUE, nrow(simulation$itempool))
    if (nrow(simulation$modules) > 0) {
      for (i in 1:nrow(simulation$modules)) {
        orphanItem[match(simulation$modules$ITEM_IDS[[i]], simulation$itempool$ITEM_ID)] = FALSE
      }
    }
    orphanItemInds = which(orphanItem)
    if (length(orphanItemInds) > 0) {
      simulation$modules = add_row(simulation$modules,
        MODULE_ID = simulation$itempool$ITEM_ID[orphanItemInds],
        ITEM_IDS = lapply(orphanItemInds, function(orphanItemInds) { return(c(simulation$itempool$ITEM_ID[orphanItemInds])) })
      )
    }
  }

  # Check if all modules are one-item modules.
  # Set a flag that can be used to enable "one-item only" performance optimizations
  if (!("oneItemModules" %in% names(simulation$control))) {
    simulation$control$oneItemModules = !any(vapply(simulation$modules$ITEM_IDS, length, as.integer(0)) > 1)
  }

  # Add 'ITEM_INDICES' to the modules, indicating the items associated with each module
  if(!"ITEM_INDICES" %in% names(simulation$modules)) {
    simulation$modules$ITEM_INDICES = lapply(simulation$modules$ITEM_IDS, function(itemIds) {
      match(itemIds, simulation$itempool$ITEM_ID)
    })
  }

  # Add 'PSG_ID' to the modules, indicating the passage associated with the items in each module
  if ("PSG_ID" %in% names(simulation$itempool)) {
    # Create an empty "PSG_ID" column, if necessary
    if (!"PSG_ID" %in% names(simulation$modules)) {
      simulation$modules$PSG_ID = as.character(NA)
    }
    noPsgIdModuleInds = which(is.na(simulation$modules$PSG_ID))
    if (any(noPsgIdModuleInds)) {
      simulation$modules$PSG_ID[noPsgIdModuleInds] = vapply(noPsgIdModuleInds, function(moduleInd) {
        # moduleInd = 1
        itemPsgIds = unique(simulation$itempool$PSG_ID[simulation$modules$ITEM_INDICES[[moduleInd]]])
        if (length(itemPsgIds) < 1) {
          return(NA)
        } else if (length(itemPsgIds) == 1) {
          return(itemPsgIds[[1]])
        } else {
          stop("module ", simulation$modules$MODULE_ID[[moduleInd]], " has too many PSG_IDs: ", paste(itemPsgIds, collapse=","))
          return(NA)
        }
      }, as.character(0))
    }
  }

  # Add 'ENEMY_CODES' to the modules, indicating modules with overlapping or conflicting items
  # This field is currently only needed by optimal item selection rules
  if((simulation$control$itemSelectionRule == "optimal") && (!"ENEMY_CODES" %in% names(simulation$modules))) {
    # initialize an empty column
    simulation$modules$ENEMY_CODES = list(character(0))
    # inspect all modules with >1 item.
    multiItemModuleInds = which(vapply(simulation$modules$ITEM_INDICES, length, as.integer(0)) > 1)
    for (moduleInd in multiItemModuleInds) {
      # moduleInd = 3
      otherModuleInds = seq_len(nrow(simulation$modules))
      otherModuleInds = otherModuleInds[otherModuleInds != moduleInd]
      for (otherModuleInd in otherModuleInds) {
        # otherModuleInd = 4
        overlapItemIds = intersect(simulation$modules$ITEM_IDS[[moduleInd]], simulation$modules$ITEM_IDS[[otherModuleInd]])
        # if any ids overlap, add the current MODULE_ID as an enemy code in both modules
        if (length(overlapItemIds) > 1) {
          simulation$modules$ENEMY_CODES[[moduleInd]] = sort(unique(c(simulation$modules$ENEMY_CODES[[moduleInd]], simulation$modules$MODULE_ID[[moduleInd]])))
          simulation$modules$ENEMY_CODES[[otherModuleInd]] = sort(unique(c(simulation$modules$ENEMY_CODES[[otherModuleInd]], simulation$modules$MODULE_ID[[moduleInd]])))
        }
      }
    }
  }

  if("constraints" %in% names(simulation)) {
    if("content" %in% names(simulation$constraints)) {
      # Add 'TYPE' to the content constraints if missing, or fill in missing values if defined
      if (!"TYPE" %in% names(simulation$constraints$content)) {
        simulation$constraints$content$TYPE = rep("H", nrow(simulation$constraints$content))
      } else {
        simulation$constraints$content$TYPE[is.na(simulation$constraints$content$TYPE)] = "H"
      }

      # Add 'CONS_INDICES' to the itempool, indicating the content constraints associated with each item
      if(!"CONS_INDICES" %in% names(simulation$itempool)) {
        simulation$itempool$CONS_INDICES = lapply(simulation$itempool$CONS_IDS, function(consIds) {
          match(consIds, simulation$constraints$content$CONS_ID)
        })
      }

      # Add 'ITEM_INDICES' to the content constraints, indicating the items associated with each content constraint
      if(!"ITEM_INDICES" %in% names(simulation$constraints$content)) {
        simulation$constraints$content$ITEM_INDICES = lapply(1:nrow(simulation$constraints$content), function(consIndex) {
          return(which(vapply(simulation$itempool$CONS_INDICES, function(itemConsIndices) {
            return(consIndex %in% itemConsIndices)
          }, as.logical(0))))
        })
      }

      # Add 'MODULE_LENGTH' to the content constraints, indicating how many items in each module associated with each content constraint
      # This field is currently only needed by optimal item selection rules
      if((simulation$control$itemSelectionRule == "optimal") && (!"MODULE_LENGTH" %in% names(simulation$constraints$content))) {
        simulation$constraints$content$MODULE_LENGTH = lapply(1:nrow(simulation$constraints$content), function(consIndex) {
          # consIndex = 1
          moduleLength = vapply(1:nrow(simulation$modules), function(moduleIndex) {
            # moduleIndex = 2
            length(intersect(simulation$constraints$content$ITEM_INDICES[[consIndex]], simulation$modules$ITEM_INDICES[[moduleIndex]]))
          }, as.integer(0))
        })
      }
    }

    if("passage" %in% names(simulation$constraints)) {
      # Add 'ITEM_INDICES' to the passage constraints, indicating the items associated with each passage
      if(!"ITEM_INDICES" %in% names(simulation$constraints$passage)) {
        simulation$constraints$passage$ITEM_INDICES = lapply(simulation$constraints$passage$PSG_ID, function(psgId) {
          # psgId = simulation$constraints$passage$PSG_ID[[1]]
          return(which(simulation$itempool$PSG_ID == psgId))
        })
      }

      # Add 'MODULE_LENGTH' to the passage constraints, indicating how many items in each module associated with each passage
      # This field is currently only needed by optimal item selection rules
      if((simulation$control$itemSelectionRule == "optimal") && (!"MODULE_LENGTH" %in% names(simulation$constraints$passage))) {
        simulation$constraints$passage$MODULE_LENGTH = lapply(1:nrow(simulation$constraints$passage), function(psgIndex) {
          # psgIndex = 1
          moduleLength = vapply(1:nrow(simulation$modules), function(moduleIndex) {
            # moduleIndex = 2
            length(intersect(simulation$constraints$passage$ITEM_INDICES[[psgIndex]], simulation$modules$ITEM_INDICES[[moduleIndex]]))
          }, as.integer(0))
        })
      }
    }
  }

  # Add 'MODULE_INDICES' to the panels, indicating the modules associated with each panel
  if (("panels" %in% names(simulation)) && !("MODULE_INDICES" %in% names(simulation$panels))) {
    simulation$panels$MODULE_INDICES = lapply(simulation$panels$MODULE_IDS, function(moduleIds) {
      match(moduleIds, simulation$modules$MODULE_ID)
    })
  }

  # Add 'PATH_INDICES' to mst$path
  if ("mst" %in% names(simulation)) {
    if (("active_path" %in% names(simulation$mst)) && !("ACTIVE_PATH_INDICES" %in% names(simulation$mst))) {
      simulation$mst$ACTIVE_PATH_INDICES = lapply(simulation$mst$active_path, function(path) {
        match(path, simulation$modules$MODULE_ID)
      })
    }
    if (("path" %in% names(simulation$mst)) && !("PATH_INDICES" %in% names(simulation$mst))) {
      simulation$mst$PATH_INDICES = lapply(simulation$mst$path, function(path) {
        match(path, simulation$modules$MODULE_ID)
      })
    }
  }

  return(simulation)
}

#' initSimulee
#'
#' @examples
#'   simulation = initSimulation(readRDS(system.file("example/passage-optimal.rds", package = "CATSimulator")))
#'   simulee = generateSimuleesByTrueTheta(-2, 10001)
#'   simuleeOut = initSimulee(simulee, simulation)
#' @import tidyverse
#' @export
initSimulee <- function (simulee, simulation) {
  suppressPackageStartupMessages(require(tidyverse))

  # Initialize the simulee... pre-allocate room for the max test length
  simuleeOut = tibble(
    MODULE_INDEX = as.integer(rep(NA, simulation$control$maxItems)),
    ITEM_INDEX = as.integer(NA),
    SCORE = as.integer(NA),
    THETA = as.double(NA),
    CSEM = as.double(NA),
    ESTIMATOR = as.character(NA),
    IEC_RELEASED = as.logical(NA),
    BLOCKED_ITEM_INDEX = list(integer(0))
  )

  # If passages are defined... create columns for PSG_ID and BLOCKED_PSG_ID
  if ("PSG_ID" %in% names(simulation$modules)) {
    simuleeOut = add_column(simuleeOut, PSG_ID = as.character(NA), .before = "MODULE_INDEX")
    simuleeOut = add_column(simuleeOut, BLOCKED_PSG_ID = list(character(0)), .before = "BLOCKED_ITEM_INDEX")
  }

  # If panels are defined, assign a PANEL_ID to the simulee
  if ("panels" %in% names(simulation)) {
    panelWeight = NULL
    if ("PANEL_WEIGHT" %in% names(simulation$panels)) {
      panelWeight = simulation$panels$PANEL_WEIGHT
    }
    simuleeOut = add_column(simuleeOut, PANEL_ID = sample(simulation$panels$PANEL_ID, size = 1, prob = panelWeight), .before = 1)
  }

  return(simuleeOut)
}

#' Run a single simulee through the given simulation.
#'
#' @param simulee A tibble containing at least one simulee, the first row will be used.
#' @param simulation An object defining the test to be run.
#' @return A tibble containing the completed simulee test output.
#' @examples
#'   options(warn=-1)
#'   simulation = initSimulation(readRDS(system.file("example/workshop2pl-small-optimal.rds", package = "CATSimulator")))
#'   simulation$control$solver = list(name = "lpsolve", external = F, mipGap = 0.0001, timeout = 1000, verbose = F)
#'   simulee = generateSimuleesByTrueTheta(-2, 10001)
#'   simuleeOut = runSimulee(simulee, simulation)
#' @import tidyverse
#' @export
runSimulee <- function(simulee, simulation) {
  suppressPackageStartupMessages(require(tidyverse))

  simuleeOut = initSimulee(simulee, simulation)

  # Select items, answer/score items, repeat... until the test is finished.
  stage = 0
  repeat {
    row = sum(!is.na(simuleeOut$ITEM_INDEX)) + 1
    stage = stage + 1

    # Set the random seed for the selection
    set.seed(simulee$SEED[[1]] + (10 * (row - 1)) + ifelse(row <= 1, 0, simuleeOut$SCORE[(row - 1)]))

    # Try to select more items
    selection = isr.select(simuleeOut, simulation, stage)
    # print(selection)
    if (is.null(selection)) {
      return(NULL)
    }
    if (is.null(selection$moduleIndex)) {
      # If no items were returned, the test is finished.
      break
    }

    # Assign the selected items to the simulee
    assignedItemInd = simulation$modules$ITEM_INDICES[[selection$moduleIndex]]
    lastRow = row + length(assignedItemInd)-1
    rowsToAssign = row:lastRow

    simuleeOut$MODULE_INDEX[rowsToAssign] = selection$moduleIndex
    simuleeOut$ITEM_INDEX[rowsToAssign] = assignedItemInd

    # Block all selected and candidate items, and all selected passage ids, so they are not eligible for future selection
    blockedItemIndex = assignedItemInd
    if (row > 1) {
      blockedItemIndex = c(blockedItemIndex, simuleeOut$BLOCKED_ITEM_INDEX[[row-1]], recursive = TRUE)
    }
    if ("candidateModuleIndices" %in% names(selection)) {
      blockedItemIndex = c(blockedItemIndex, simulation$modules$ITEM_INDICES[selection$candidateModuleIndices], recursive = TRUE)
    }
    simuleeOut$BLOCKED_ITEM_INDEX[[lastRow]] = sort(unique(blockedItemIndex))

    if (("PSG_ID" %in% names(simulation$modules)) && ("psgId" %in% names(selection))) {
      simuleeOut$PSG_ID[rowsToAssign] = selection$psgId

      blockedPsgIds = selection$psgId
      if (row > 1) {
        blockedPsgIds = c(blockedPsgIds, simuleeOut$BLOCKED_PSG_ID[[row-1]])
      }
      simuleeOut$BLOCKED_PSG_ID[[lastRow]] = sort(unique(blockedPsgIds))
    }

    # Track if and where IEC blocked items were released. Do not clear the blockedItemIndices, we keep collecting but ignore them.
    iecReleased = ("iecReleased" %in% names(selection)) && selection$iecReleased
    if (row > 1) {
      iecReleased = iecReleased || simuleeOut$IEC_RELEASED[[row-1]]
    }
    simuleeOut$IEC_RELEASED[[lastRow]] = iecReleased

    # Answer/score each item that was selected.
    simuleeOut$SCORE[rowsToAssign] = generateScores(simulation$item_pars[simuleeOut$ITEM_INDEX[rowsToAssign],,drop=FALSE], simulee$TRUE_THETA[[1]])

    # Calculate theta/csem based on all assigned items/scores
    rowsToEstimate = 1:lastRow
    item_parsToEstimate = simulation$item_pars[simuleeOut$ITEM_INDEX[rowsToEstimate],,drop=FALSE]
    scoresToEstimate = simuleeOut$SCORE[rowsToEstimate]
    abilityEstimate = estimateAbility(item_parsToEstimate, scoresToEstimate, simulation)
    simuleeOut$THETA[lastRow] = abilityEstimate$theta
    simuleeOut$CSEM[lastRow] = abilityEstimate$csem
    simuleeOut$ESTIMATOR[lastRow] = abilityEstimate$abilityEstimator
  }

  return (simuleeOut)
}

#' Generate simulees and run them through the given simulation.
#'
#' @param simulation An object defining the test to be run.
#' @param simulees A tibble of simulees defining id, true theta, and random number seed
#' @param progressCallback A callback function that takes a single parameter, the number of completed simulees.
#' @param show_progress A boolean indicating whether to write progress events to the console
#' @return A tibble of simulee test output.
#' @examples
#'   options(warn=-1)
#'   simulation = initSimulation(readRDS(system.file("example/workshop2pl-small-optimal.rds", package = "CATSimulator")))
#'   simulation$control$solver = list(name = "lpsolve", external = F, mipGap = 0.0001, timeout = 1000, verbose = F)
#'   simulees = generateSimuleesByTrueTheta(-2:2, 10001:10005)
#'   simulationOut = runSimulation(simulation, simulees)
#' @import tidyverse
#' @export
runSimulation <- function(simulation, simulees, progressCallback = NULL, showProgress = F) {
  if (is.null(simulation)) {
    return(NULL)
  }
  suppressPackageStartupMessages(require(tidyverse))

  simulationOut = tibble(
    SIM_ID = character(),
    SEED = integer(),
    TRUE_THETA = double(),
    MODULE_ID = character(),
    ITEM_ID = character(),
    SCORE = integer(),
    THETA = double(),
    CSEM = double(),
    ESTIMATOR = character()
  )
  if ("PSG_ID" %in% names(simulation$modules)) {
    simulationOut = add_column(simulationOut, .after = "MODULE_ID",
                               PSG_ID = character()
    )
  }

  if (showProgress) {
    cat("Simulation start at ", format(Sys.time(), "%X"), "\n")
  }

  for (simuleeInd in 1:nrow(simulees)) {
    # Run the simulee
    simuleeOut = runSimulee(simulees[simuleeInd,], simulation)

    # Append the completed rows of the simuleeOut to the combined simulationOut
    rowsToAppend = 1:sum(!is.na(simuleeOut$ITEM_INDEX))
    simulationOut = add_row(simulationOut,
                            SIM_ID = simulees$ID[[simuleeInd]],
                            SEED = simulees$SEED[[simuleeInd]],
                            TRUE_THETA = simulees$TRUE_THETA[[simuleeInd]],
                            MODULE_ID = simulation$modules$MODULE_ID[simuleeOut$MODULE_INDEX[rowsToAppend]],
                            ITEM_ID = simulation$itempool$ITEM_ID[simuleeOut$ITEM_INDEX[rowsToAppend]],
                            SCORE = simuleeOut$SCORE[rowsToAppend],
                            THETA = simuleeOut$THETA[rowsToAppend],
                            CSEM = simuleeOut$CSEM[rowsToAppend],
                            ESTIMATOR = simuleeOut$ESTIMATOR[rowsToAppend]
    )
    if ("PSG_ID" %in% names(simulation$modules)) {
      simulationOut$PSG_ID[simulationOut$SIM_ID == simulees$ID[[simuleeInd]]] = simuleeOut$PSG_ID[rowsToAppend]
    }

    # If there is a progress callback, report progress; otherwise, printout progress on console
    if (!is.null(progressCallback)) {
      progressCallback(simuleeInd)
    }
    if (showProgress && (simuleeInd %% 10 == 1)) {
      cat(".")
    }
  }

  # Join the item pars immediately after the ITEM_ID
  simulationOut <- simulationOut %>%
    left_join(simulation$itempool %>% select(!any_of(c("CONS_IDS", "CONS_INDICES"))), by = "ITEM_ID") %>%
    select("SIM_ID", "SEED", "TRUE_THETA", "MODULE_ID", "ITEM_ID", any_of(c("MODEL", "NC")), starts_with("PAR_"), everything())

  if (showProgress) {
    cat("Simulation end at ", format(Sys.time(), "%X"), "\n")
  }
  return(simulationOut)
}

#' Generate simulees and run them in parallel through the given simulation.
#'
#' The caller is responsible for setup, management and cleanup of the future::plan().
#'     future::plan(future::multisession, workers = future::availableCores(omit = 1))
#'
#' @param simulation An object defining the test to be run.
#' @param simulees A tibble of simulees defining id, true theta, and random number seed
#' @param progressCallback A callback function that takes a single parameter, the number of completed simulees.
#' @param show_progress A boolean indicating whether to write progress events to the console
#' @return A tibble of simulee test output.
#' @examples
#'   options(warn=-1)
#'   simulation = initSimulation(readRDS(system.file("example/workshop2pl-small-optimal.rds", package = "CATSimulator")))
#'   simulation$control$solver = list(name = "lpsolve", external = F, mipGap = 0.0001, timeout = 1000, verbose = F)
#'   simulees = generateSimuleesByTrueTheta(-2:2, 10001:10005)
#'   simulationOut = runParallelSimulation(simulation, simulees)
#' @import tidyverse
#' @export
runParallelSimulation <- function(simulation, simulees, progressCallback = NULL, showProgress = F) {
  if (is.null(simulation)) {
    return(NULL)
  }
  suppressPackageStartupMessages(require(tidyverse))
  if (!require("furrr", quietly = TRUE)) {
    stop("Package \"furrr\" needed for this function to work. Please install it or use 'runSimulation' instead.", call. = FALSE)
  }

  simulationOut = tibble(
    SIM_ID = character(),
    SEED = integer(),
    TRUE_THETA = double(),
    MODULE_ID = character(),
    ITEM_ID = character(),
    SCORE = integer(),
    THETA = double(),
    CSEM = double(),
    ESTIMATOR = character()
  )
  if ("PSG_ID" %in% names(simulation$modules)) {
    simulationOut = add_column(simulationOut, .after = "MODULE_ID", PSG_ID = character())
  }

  if (showProgress) {
    cat("Simulation start at ", format(Sys.time(), "%X"), "\n")
  }

  simulationOut <- simulees %>%
    furrr::future_pmap(function(...) {
      # each future is its own environment with its own imports
      suppressPackageStartupMessages(require(tidyverse))

      # Run the simulee
      simulee = tibble_row(...)
      simuleeOut = runSimulee(simulee, simulation)

      # Trim unused rows
      simuleeOut = simuleeOut[!is.na(simuleeOut$ITEM_INDEX),]

      # Transform to simulationOut columns
      simuleeOut = simuleeOut %>%
        mutate(MODULE_ID = simulation$modules$MODULE_ID[MODULE_INDEX]) %>%
        mutate(ITEM_ID = simulation$itempool$ITEM_ID[ITEM_INDEX]) %>%
        select(matches("PSG_ID"), MODULE_ID, ITEM_ID, SCORE, THETA, CSEM, ESTIMATOR)

      # If there is a progress callback, report progress; otherwise, printout progress on console
      simuleeInd = which(simulees$ID == simulee$ID)
      if (!is.null(progressCallback)) {
        progressCallback(simuleeInd)
      }
      if (showProgress && (simuleeInd %% 10 == 1)) {
        cat(".")
      }

      return(bind_cols(simulee %>% rename(SIM_ID = ID), simuleeOut))
    }, .options = furrr_options(
      # globals = list(future.globals.maxSize = Inf),
      seed = lapply(simulees$SEED, function(x) {
        # see https://rdrr.io/github/DavisVaughan/furrr/src/R/seed.R
        # function as_lecyer_cmrg_seed_from_integer
        RNGkind("L'Ecuyer-CMRG")
        set.seed(x)
        return(.Random.seed)
      })
    )) %>%
    bind_rows()

  # Join the item pars immediately after the ITEM_ID
  simulationOut <- simulationOut %>%
    left_join(simulation$itempool %>% select(!c("CONS_IDS", "CONS_INDICES")), by = "ITEM_ID") %>%
    select("SIM_ID", "SEED", "TRUE_THETA", "MODULE_ID", "ITEM_ID", any_of(c("MODEL", "NC")), starts_with("PAR_"), everything())

  if (showProgress) {
    cat("Simulation end at ", format(Sys.time(), "%X"), "\n")
  }

  return(simulationOut)
}