library(RcppAlgos)
library(fmsb)
library(matrixStats)

N <- 50

# Generate all combinations
cm <- comboGeneral(0:N, 3, repetition = FALSE)

# Pre-allocate
n_vec <- rep(N, 3)
success <- cm
failure <- matrix(n_vec, nrow = nrow(cm), ncol = 3, byrow = TRUE) - cm

# Modified fast Fisher function - correct implementation
fast_fisher <- function(success, failure) {
# Initialize p-value matrix
pvals <- matrix(NA_real_, nrow = nrow(success), ncol = 3)

# Perform pairwise comparisons
for (i in 1:nrow(success)) {
# 1 vs 2
pvals[i,1] <- fisher.test(matrix(c(success[i,1], failure[i,1],
success[i,2], failure[i,2]), nrow = 2))$p.value

# 1 vs 3
pvals[i,2] <- fisher.test(matrix(c(success[i,1], failure[i,1],
success[i,3], failure[i,3]), nrow = 2))$p.value

# 2 vs 3
pvals[i,3] <- fisher.test(matrix(c(success[i,2], failure[i,2],
success[i,3], failure[i,3]), nrow = 2))$p.value
}

# Bonferroni adjustment
pmin(pvals * 3, 1) # Cap at 1 after adjustment
}

# Run with timing
system.time({
# Overall Fisher tests
overall_p <- apply(cbind(success, failure), 1, function(x) {
fisher.test(matrix(x, nrow = 2))$p.value
})

# Pairwise Fisher tests
pairwise_p <- fast_fisher(success, failure)
min_pairwise_p <- rowMins(pairwise_p, na.rm = TRUE)

# Filter condition
keep <- overall_p > 0.05 & min_pairwise_p < 0.05
result <- cm[keep, ]
})

# Print first few results
head(result)
nrow(result) # Number of qualifying combinations