library(parallel)
library(data.table)

# Define parameters
Part <- 40
REP <- 1000
thresholds <- c(0, 0.025, 0.05, 0.075, 0.1)

# Define functions
pvaluefromreplace <- function(M, C, D, ql, qh, replace = TRUE) {
  if (ql) QL <- quantile(D[[M]], ql) else QL <- min(D[[M]]) - 1
  if (qh) QH <- quantile(D[[M]], 1 - qh) else QH <- max(D[[M]]) + 1
  
  if (!replace) {
    DD <- D[D[[M]] > QL & D[[M]] < QH]
  } else {
    DD <- copy(D)
    DD[[M]][DD[[M]] < QL] <- QL
    DD[[M]][DD[[M]] > QH] <- QH
  }
  
  a <- summary(aov(DD[[M]] ~ DD[[C]]))
  return(a[[1]][["Pr(>F)"]][1])
}

pvaluefrom <- function(M, C, D, R) {
  DD <- D[R]
  a <- summary(aov(DD[[M]] ~ DD[[C]]))
  return(a[[1]][["Pr(>F)"]][1])
}

# Function to summarize p-values
summarize_pvalues <- function(data, restriction) {
  filtered_data <- data[eval(substitute(restriction), data), ]
  
  # Calculate the summary statistics
  summary <- data.frame(
    pvalue_below_05 = sum(filtered_data$pvalue < 0.05, na.rm = TRUE),
    pvalue_above_05 = sum(filtered_data$pvalue >= 0.05, na.rm = TRUE),
    quantile_pvalue_05 = quantile(filtered_data$pvalue, 0.05, na.rm = TRUE)
  )
  summary$FER = summary$pvalue_below_05 / (summary$pvalue_below_05+summary$pvalue_above_05)
  
  return(summary)
}

run_simulation <- function(x) {
  # Simulate some data
  X <- data.table(
    Cond = rep(1:2, each = Part / 2),
    RT = rnorm(Part),
    EXCL = rnorm(Part)
  )
  
  # Create all combinations of parameters
  sim_results <- expand.grid(
    simulation = x,
    replacement = c(TRUE, FALSE),
    differentmeasure = c(FALSE, TRUE),
    lowexclusion = thresholds,
    highexclusion = thresholds
  )
  
  # Filter combinations:
  sim_results <- sim_results[!(sim_results$replacement == TRUE & sim_results$differentmeasure == TRUE), ]
  
  # Add pvalue calculations
  sim_results$pvalue <- mapply(
    function(replace, diff, ql, qh) {
      if (replace & !diff) {
        pvaluefromreplace("RT", "Cond", X, ql, qh, replace = TRUE)
      } else if (!replace & !diff) {
        pvaluefromreplace("RT", "Cond", X, ql, qh, replace = FALSE)
      } else {
        pvaluefrom("RT", "Cond", X, X$EXCL > quantile(X$EXCL, ql) & X$EXCL < quantile(X$EXCL, 1 - qh))
      }
    },
    sim_results$replacement,
    sim_results$differentmeasure,
    sim_results$lowexclusion,
    sim_results$highexclusion
  )
  
  return(as.data.table(sim_results))
}

# Use mclapply
RESULTS <- mclapply(1:REP, run_simulation, mc.cores = detectCores() - 1)

# Combine all simulation results
RESULTS <- rbindlist(RESULTS, fill = TRUE)

# Single threshold
single_summary <- summarize_pvalues(RESULTS, lowexclusion == 0 & highexclusion == 0 & replacement == FALSE & differentmeasure==FALSE)

# Multiple thresholds, without replacement
multi_no_replace <- RESULTS[
  (replacement == FALSE & differentmeasure == FALSE) | (lowexclusion == 0 & highexclusion == 0), 
]
multi_no_replace <- multi_no_replace[, .(pvalue = min(pvalue)), by = simulation]
multi_no_replace_summary <- summarize_pvalues(multi_no_replace, TRUE)

# Multiple thresholds, with replacement
multi_with_replace <- RESULTS[replacement == TRUE & differentmeasure == FALSE, ]
multi_with_replace <- multi_with_replace[, .(pvalue = min(pvalue)), by = simulation]
multi_with_replace_summary <- summarize_pvalues(multi_with_replace, TRUE)

# Multiple thresholds, different measure (without replacement)
multi_diff_no_replace <- RESULTS[
  (replacement == FALSE & differentmeasure == TRUE) | (lowexclusion == 0 & highexclusion == 0), 
]
multi_diff_no_replace <- multi_diff_no_replace[, .(pvalue = min(pvalue)), by = simulation]
multi_diff_no_replace_summary <- summarize_pvalues(multi_diff_no_replace, TRUE)
# Combine all summaries
summary_table <- rbind(
  cbind(type = "No threshold", single_summary),
  cbind(type = "Multiple thresholds (no replacement)", multi_no_replace_summary),
  cbind(type = "Multiple thresholds (with replacement)", multi_with_replace_summary),
  cbind(type = "Multiple thresholds (different measure, no replacement)", multi_diff_no_replace_summary)
)

# Display the final table
print(summary_table)