Performance Benchmarks

# Create test data
n_ct <- 10
n_genes <- 100
test_data <- matrix(runif(n_ct * n_genes), nrow = n_ct)

# Compare implementations
n_iter <- 50
times_r <- numeric(n_iter)
times_cpp <- numeric(n_iter)

for (i in 1:n_iter) {
  t1 <- system.time({
    compute_correlation(test_data, use_cpp = FALSE)
  })
  times_r[i] <- t1["elapsed"]
  
  t2 <- system.time({
    compute_correlation(test_data, use_cpp = TRUE)
  })
  times_cpp[i] <- t2["elapsed"]
}

df_bench <- data.frame(
  Time = c(times_r, times_cpp) * 1000,
  Implementation = rep(c("R (vectorized)", "C++ (RcppArmadillo)"), each = n_iter)
)

ggplot(df_bench, aes(x = Implementation, y = Time, fill = Implementation)) +
  geom_boxplot(alpha = 0.7) +
  scale_fill_manual(values = c("#e74c3c", "#3498db")) +
  labs(
    title = "R vs C++ Implementation Performance",
    subtitle = paste("Computing correlation for", n_ct, "× ", n_genes, "matrix"),
    x = "",
    y = "Time (milliseconds)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

gene_sizes <- c(100, 200, 500, 1000, 2000)
n_ct <- 8
results <- list()

for (n_genes in gene_sizes) {
  test_data <- matrix(runif(n_ct * n_genes), nrow = n_ct)
  
  # Time correlation computation
  t <- system.time({
    for (j in 1:20) {
      compute_correlation(test_data, use_cpp = TRUE)
    }
  })
  
  results[[length(results) + 1]] <- data.frame(
    n_genes = n_genes,
    time = t["elapsed"] / 20 * 1000
  )
}

df_scale_genes <- do.call(rbind, results)

ggplot(df_scale_genes, aes(x = n_genes, y = time)) +
  geom_line(color = "#3498db", linewidth = 1) +
  geom_point(color = "#3498db", size = 3) +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    title = "Performance Scaling with Gene Count",
    subtitle = paste(n_ct, "cell types"),
    x = "Number of Genes (log scale)",
    y = "Time per iteration (ms, log scale)"
  ) +
  theme_minimal(base_size = 12)

ct_sizes <- c(3, 5, 10, 20, 30)
n_genes <- 500
results_ct <- list()

for (n_ct in ct_sizes) {
  test_data <- matrix(runif(n_ct * n_genes), nrow = n_ct)
  
  t <- system.time({
    for (j in 1:20) {
      compute_correlation(test_data, use_cpp = TRUE)
    }
  })
  
  results_ct[[length(results_ct) + 1]] <- data.frame(
    n_ct = n_ct,
    time = t["elapsed"] / 20 * 1000
  )
}

df_scale_ct <- do.call(rbind, results_ct)

ggplot(df_scale_ct, aes(x = n_ct, y = time)) +
  geom_line(color = "#e74c3c", linewidth = 1) +
  geom_point(color = "#e74c3c", size = 3) +
  labs(
    title = "Performance Scaling with Cell Type Count",
    subtitle = paste(n_genes, "genes"),
    x = "Number of Cell Types",
    y = "Time per iteration (ms)"
  ) +
  theme_minimal(base_size = 12)

# Setup
n_ct <- 8
n_genes <- 1000
reference <- matrix(abs(rnorm(n_ct * n_genes, 2)), nrow = n_ct)
rownames(reference) <- paste0("CT", 1:n_ct)
colnames(reference) <- paste0("Gene", 1:n_genes)

# Benchmark components
timings <- list()

# Initialization
t_init <- system.time({
  dw <- darwin(reference)
})
timings$Initialization <- t_init["elapsed"]

# Optimization (short run for benchmark)
t_opt <- system.time({
  dw$optimize(ngen = 20, pop_size = 50, verbose = FALSE, parallel = FALSE)
})
timings$Optimization <- t_opt["elapsed"]

# Selection
t_sel <- system.time({
  dw$select(weights = c(-1, 1))
})
timings$Selection <- t_sel["elapsed"]

# Create timing data frame
df_timing <- data.frame(
  Step = names(timings),
  Time = unlist(timings)
)
df_timing$Percent <- df_timing$Time / sum(df_timing$Time) * 100

ggplot(df_timing, aes(x = reorder(Step, Time), y = Time, fill = Step)) +
  geom_bar(stat = "identity", alpha = 0.8) +
  geom_text(aes(label = paste0(round(Percent, 1), "%")), hjust = -0.1) +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Workflow Time Breakdown",
    subtitle = paste(n_genes, "genes,", n_ct, "cell types, 20 generations"),
    x = "",
    y = "Time (seconds)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  expand_limits(y = max(df_timing$Time) * 1.2)

gen_sizes <- c(10, 25, 50, 100)
results_gen <- list()

for (ngen in gen_sizes) {
  t <- system.time({
    dw <- darwin(reference)
    dw$optimize(ngen = ngen, pop_size = 50, verbose = FALSE, parallel = FALSE)
  })
  
  results_gen[[length(results_gen) + 1]] <- data.frame(
    generations = ngen,
    time = t["elapsed"]
  )
}

df_scale_gen <- do.call(rbind, results_gen)

# Fit linear model
lm_fit <- lm(time ~ generations, data = df_scale_gen)

ggplot(df_scale_gen, aes(x = generations, y = time)) +
  geom_line(color = "#3498db", linewidth = 1) +
  geom_point(color = "#3498db", size = 3) +
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed", color = "gray50") +
  labs(
    title = "Scaling with Number of Generations",
    subtitle = paste("~", round(coef(lm_fit)[2], 3), "seconds per generation"),
    x = "Number of Generations",
    y = "Total Time (seconds)"
  ) +
  theme_minimal(base_size = 12)

pop_sizes <- c(30, 50, 80, 100, 150)
results_pop <- list()

for (pop in pop_sizes) {
  t <- system.time({
    dw <- darwin(reference)
    dw$optimize(ngen = 30, pop_size = pop, verbose = FALSE, parallel = FALSE)
  })
  
  results_pop[[length(results_pop) + 1]] <- data.frame(
    pop_size = pop,
    time = t["elapsed"]
  )
}

df_scale_pop <- do.call(rbind, results_pop)

ggplot(df_scale_pop, aes(x = pop_size, y = time)) +
  geom_line(color = "#e74c3c", linewidth = 1) +
  geom_point(color = "#e74c3c", size = 3) +
  labs(
    title = "Scaling with Population Size",
    subtitle = "30 generations",
    x = "Population Size",
    y = "Total Time (seconds)"
  ) +
  theme_minimal(base_size = 12)