Run Time Recording

# Load relevant libraries
library(SurprisalAnalysis)
library(ggplot2)
library(peakRAM)

Please uncomment the below and run

#setting seed to ensure reproducibility
set.seed(123)

#params
n_genes       <- 1000
sample_counts <- c(5, 10, 20, 50, 100)
n_reps        <- 3


#' Function records run time and RAM peak usage
#'
#' @param n_genes number of genes
#' @param n_samples
#'
#' @return run time and RAM peak usage
#' @export
time_ram_single_run <- function(n_genes, n_samples) {
  elapsed <- NA_real_
  pr <- peakRAM({

    expr_matrix <- matrix(
      rlnorm(n = n_genes * n_samples, meanlog = 2, sdlog = 1),
      nrow = n_genes,
      ncol = n_samples
    )
    expr_df <- as.data.frame(expr_matrix)
    rownames(expr_df) <- paste0("Gene", seq_len(n_genes))
    colnames(expr_df) <- paste0("Sample", seq_len(n_samples))

    t0 <- Sys.time()
    res <- surprisal_analysis(expr_df)
    t1 <- Sys.time()
    elapsed <- as.numeric(difftime(t1, t0, units = "secs"))
  })
  list(
    time_sec = elapsed,
    peak_ram_mib = pr$Peak_RAM_Used_MiB[1]
  )
}


#store info

results_df <- data.frame(
  n_samples               = integer(length(sample_counts)),
  mean_time               = numeric(length(sample_counts)),
  sd_time                 = numeric(length(sample_counts)),
  ci_lower                = numeric(length(sample_counts)),
  ci_upper                = numeric(length(sample_counts)),
  mean_peak_ram_mib       = numeric(length(sample_counts)),
  sd_peak_ram_mib         = numeric(length(sample_counts)),
  ci_lower_peak_ram_mib   = numeric(length(sample_counts)),
  ci_upper_peak_ram_mib   = numeric(length(sample_counts)),
  stringsAsFactors = FALSE
)


ci_halfwidth_fun <- function(x, n) {
  sdev <- sd(x)
  se   <- sdev / sqrt(n)
  qt(0.975, df = n - 1) * se
}

for (i in seq_along(sample_counts)) {
  s <- sample_counts[i]
  rep_times <- numeric(n_reps)
  rep_ram   <- numeric(n_reps)

  for (r in seq_len(n_reps)) {
    gc()
    tr <- time_ram_single_run(n_genes = n_genes, n_samples = s)
    rep_times[r] <- tr$time_sec
    rep_ram[r]   <- tr$peak_ram_mib
  }

  mu_time <- mean(rep_times)
  sd_time <- sd(rep_times)
  hw_time <- ci_halfwidth_fun(rep_times, n_reps)

  mu_ram  <- mean(rep_ram)
  sd_ram  <- sd(rep_ram)
  hw_ram  <- ci_halfwidth_fun(rep_ram, n_reps)

  results_df[i, ] <- list(
    s,
    mu_time, sd_time, mu_time - hw_time, mu_time + hw_time,
    mu_ram,  sd_ram,  mu_ram  - hw_ram,  mu_ram  + hw_ram
  )
  print(i)
  message(sprintf(
    "n_samples = %3d : time = %.4f s (sd=%.4f, 95%% CI=[%.4f, %.4f]);  peak RAM = %.2f MiB (sd=%.2f, 95%% CI=[%.2f, %.2f])",
    s, mu_time, sd_time, mu_time - hw_time, mu_time + hw_time,
    mu_ram, sd_ram, mu_ram - hw_ram, mu_ram + hw_ram
  ))
}

## [1] 1

## n_samples =   5 : time = 0.0020 s (sd=0.0000, 95% CI=[0.0020, 0.0021]);  peak RAM = 0.60 MiB (sd=0.00, 95% CI=[0.60, 0.60])

## [1] 2

## n_samples =  10 : time = 0.0047 s (sd=0.0001, 95% CI=[0.0045, 0.0050]);  peak RAM = 2.40 MiB (sd=0.00, 95% CI=[2.40, 2.40])

## [1] 3

## n_samples =  20 : time = 0.0180 s (sd=0.0004, 95% CI=[0.0171, 0.0189]);  peak RAM = 11.80 MiB (sd=0.00, 95% CI=[11.80, 11.80])

## [1] 4

## n_samples =  50 : time = 0.1922 s (sd=0.0171, 95% CI=[0.1497, 0.2347]);  peak RAM = 199.80 MiB (sd=0.00, 95% CI=[199.80, 199.80])

## [1] 5

## n_samples = 100 : time = 1.4288 s (sd=0.0275, 95% CI=[1.3606, 1.4970]);  peak RAM = 382.00 MiB (sd=0.00, 95% CI=[382.00, 382.00])

print(results_df)

##   n_samples   mean_time      sd_time    ci_lower    ci_upper mean_peak_ram_mib
## 1         5 0.002014716 0.0000240275 0.001955029 0.002074404               0.6
## 2        10 0.004723867 0.0001101644 0.004450203 0.004997530               2.4
## 3        20 0.018030882 0.0003588299 0.017139499 0.018922265              11.8
## 4        50 0.192215761 0.0171117200 0.149707892 0.234723630             199.8
## 5       100 1.428794066 0.0274699237 1.360554993 1.497033140             382.0
##   sd_peak_ram_mib ci_lower_peak_ram_mib ci_upper_peak_ram_mib
## 1               0                   0.6                   0.6
## 2               0                   2.4                   2.4
## 3               0                  11.8                  11.8
## 4               0                 199.8                 199.8
## 5               0                 382.0                 382.0

model_cubic <- lm(mean_time ~ poly(n_samples, 3), data = results_df)
summary(model_cubic)

## 
## Call:
## lm(formula = mean_time ~ poly(n_samples, 3), data = results_df)
## 
## Residuals:
##          1          2          3          4          5 
##  4.586e-05 -8.169e-05  4.085e-05 -5.446e-06  4.300e-07 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.292e-01  4.577e-05  7191.3 8.85e-05 ***
## poly(n_samples, 3)1 1.170e+00  1.023e-04 11436.1 5.57e-05 ***
## poly(n_samples, 3)2 4.068e-01  1.023e-04  3974.3  0.00016 ***
## poly(n_samples, 3)3 3.959e-02  1.023e-04   386.8  0.00165 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0001023 on 1 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 4.891e+07 on 3 and 1 DF,  p-value: 0.0001051

results_df$pred_cubic <- predict(model_cubic, newdata = results_df)

ggplot(results_df, aes(x = n_samples, y = mean_time)) +
  geom_ribbon(aes(ymin = mean_time - sd_time, ymax = mean_time + sd_time),
              fill = "steelblue", alpha = 0.2) +
  geom_point(color = "steelblue4", size = 2, shape=8, stroke=1.5) +
  geom_line(color = "steelblue4") +
  stat_smooth(method = "lm",
              formula = y ~ poly(x, 3),
              se = FALSE,
              color = "firebrick",
              linetype = "dashed",
              size = 1) +
  labs(
    x = "Number of samples",
    y = "Average elapsed time (seconds)",
    title = sprintf("SurprisalAnalysis Runtime vs. Number of Samples\n(gene count = %d; mean ± 95%% CI over %d runs)", n_genes, n_reps), tag="A"
  ) +
  theme_minimal(base_size = 12)

is_bytes <- max(results_df$mean_peak_ram_mib, na.rm = TRUE) > 1e5
scale_factor <- if (is_bytes) 1/(1024^2) else 1

use_ci <- all(c("ci_lower_peak_ram_mib", "ci_upper_peak_ram_mib") %in% names(results_df))

df_plot <- within(results_df, {
  mean_ram_plot  <- mean_peak_ram_mib * scale_factor
  lower_ram_plot <- (if (use_ci) ci_lower_peak_ram_mib else mean_peak_ram_mib - sd_peak_ram_mib) * scale_factor
  upper_ram_plot <- (if (use_ci) ci_upper_peak_ram_mib else mean_peak_ram_mib + sd_peak_ram_mib) * scale_factor
})


ggplot(df_plot, aes(x = n_samples, y = mean_ram_plot)) +
  geom_ribbon(aes(ymin = lower_ram_plot, ymax = upper_ram_plot),
              fill = "steelblue", alpha = 0.2) +
  geom_point(color = "steelblue4", size = 2, shape=8, stroke = 1.5) +
  geom_line(color = "steelblue4") +
  scale_x_continuous(breaks = sample_counts) +
  labs(
    title = sprintf("SurprisalAnalysis Peak RAM vs. Number of Samples\n(gene count = %d; mean ± 95%% CI over %d runs)", n_genes, n_reps),
    x = "Number of Samples",
    y = "Peak RAM (MiB)", tag="B"
  ) +
  theme_minimal(base_size = 12)

#write.csv(results_df, "~/Downloads/runtime_altering_samples_try_2.csv")

####Altering the number of genes


#params
n_samples_fixed <- 50
gene_counts     <- c(100, 500, 1000, 2000, 5000)
n_reps          <- 3


results_genes_df <- data.frame(
  n_genes                 = integer(length(gene_counts)),
  mean_time               = numeric(length(gene_counts)),
  sd_time                 = numeric(length(gene_counts)),
  ci_lower                = numeric(length(gene_counts)),
  ci_upper                = numeric(length(gene_counts)),
  mean_peak_ram_mib       = numeric(length(gene_counts)),
  sd_peak_ram_mib         = numeric(length(gene_counts)),
  ci_lower_peak_ram_mib   = numeric(length(gene_counts)),
  ci_upper_peak_ram_mib   = numeric(length(gene_counts)),
  stringsAsFactors = FALSE
)

ci_halfwidth_fun <- function(x, n) {
  sdev <- sd(x); se <- sdev / sqrt(n); qt(0.975, df = n - 1) * se
}


for (i in seq_along(gene_counts)) {
  g <- gene_counts[i]
  rep_times <- numeric(n_reps)
  rep_ram   <- numeric(n_reps)

  for (r in seq_len(n_reps)) {
    gc()
    tr <- time_ram_single_run(n_genes = g, n_samples = n_samples_fixed)
    rep_times[r] <- tr$time_sec
    rep_ram[r]   <- tr$peak_ram_mib
  }

  mu_time <- mean(rep_times); sd_time <- sd(rep_times); hw_time <- ci_halfwidth_fun(rep_times, n_reps)
  mu_ram  <- mean(rep_ram);   sd_ram  <- sd(rep_ram);  hw_ram  <- ci_halfwidth_fun(rep_ram,   n_reps)

  results_genes_df[i, ] <- list(
    g,
    mu_time, sd_time, mu_time - hw_time, mu_time + hw_time,
    mu_ram,  sd_ram,  mu_ram  - hw_ram,  mu_ram  + hw_ram
  )

  message(sprintf(
    "n_genes = %5d : time = %.4f s (sd=%.4f, 95%% CI=[%.4f, %.4f]);  peak RAM = %.2f MiB (sd=%.2f, 95%% CI=[%.2f, %.2f])",
    g, mu_time, sd_time, mu_time - hw_time, mu_time + hw_time,
    mu_ram, sd_ram, mu_ram - hw_ram, mu_ram + hw_ram
  ))
}

## n_genes =   100 : time = 0.1757 s (sd=0.0039, 95% CI=[0.1660, 0.1854]);  peak RAM = 163.90 MiB (sd=0.00, 95% CI=[163.90, 163.90])

## n_genes =   500 : time = 0.1771 s (sd=0.0018, 95% CI=[0.1726, 0.1816]);  peak RAM = 179.90 MiB (sd=0.00, 95% CI=[179.90, 179.90])

## n_genes =  1000 : time = 0.1808 s (sd=0.0028, 95% CI=[0.1739, 0.1878]);  peak RAM = 199.80 MiB (sd=0.00, 95% CI=[199.80, 199.80])

## n_genes =  2000 : time = 0.1868 s (sd=0.0011, 95% CI=[0.1839, 0.1896]);  peak RAM = 239.70 MiB (sd=0.00, 95% CI=[239.70, 239.70])

## n_genes =  5000 : time = 0.2115 s (sd=0.0033, 95% CI=[0.2032, 0.2198]);  peak RAM = 359.50 MiB (sd=0.00, 95% CI=[359.50, 359.50])

print(results_genes_df)

##   n_genes mean_time     sd_time  ci_lower  ci_upper mean_peak_ram_mib
## 1     100 0.1756926 0.003917552 0.1659609 0.1854244             163.9
## 2     500 0.1770931 0.001813520 0.1725881 0.1815981             179.9
## 3    1000 0.1808388 0.002787599 0.1739140 0.1877636             199.8
## 4    2000 0.1867924 0.001148495 0.1839394 0.1896454             239.7
## 5    5000 0.2115283 0.003345237 0.2032183 0.2198383             359.5
##   sd_peak_ram_mib ci_lower_peak_ram_mib ci_upper_peak_ram_mib
## 1               0                 163.9                 163.9
## 2               0                 179.9                 179.9
## 3               0                 199.8                 199.8
## 4               0                 239.7                 239.7
## 5               0                 359.5                 359.5

ggplot(results_genes_df, aes(x = n_genes, y = mean_time)) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
              fill = "steelblue", alpha = 0.2) +
  geom_point(color = "steelblue4", size = 2, shape = 8, stroke = 1.5) +
  geom_line(color = "steelblue4") +
  geom_smooth(method = "lm", se = FALSE, color = "firebrick", size = 1) +
  scale_x_continuous(breaks = gene_counts) +
  labs(
    x = "Number of genes",
    y = "Average elapsed time (seconds)",
    title = sprintf(
      "SurprisalAnalysis Runtime vs. Number of Genes\n(sample count = %d; mean ± 95%% CI over %d runs)",
      n_samples_fixed, n_reps
    ),
    tag = "C"
  ) +
  theme_minimal(base_size = 12)

## `geom_smooth()` using formula = 'y ~ x'

is_bytes_g <- max(results_genes_df$mean_peak_ram_mib, na.rm = TRUE) > 1e5
scale_factor_g <- if (is_bytes_g) 1/(1024^2) else 1
use_ci_g <- all(c("ci_lower_peak_ram_mib", "ci_upper_peak_ram_mib") %in% names(results_genes_df))

df_plot_genes <- within(results_genes_df, {
  mean_ram_plot_g  <- mean_peak_ram_mib * scale_factor_g
  lower_ram_plot_g <- (if (use_ci_g) ci_lower_peak_ram_mib else mean_peak_ram_mib - sd_peak_ram_mib) * scale_factor_g
  upper_ram_plot_g <- (if (use_ci_g) ci_upper_peak_ram_mib else mean_peak_ram_mib + sd_peak_ram_mib) * scale_factor_g
})


ggplot(df_plot_genes, aes(x = n_genes, y = mean_ram_plot_g)) +
  geom_ribbon(aes(ymin = lower_ram_plot_g, ymax = upper_ram_plot_g),
              fill = "steelblue", alpha = 0.2) +
  geom_point(color = "steelblue4", size = 2, shape = 8, stroke = 1.5) +
  geom_line(color = "steelblue4") +
  scale_x_continuous(breaks = gene_counts) +
  labs(
    title = sprintf("SurprisalAnalysis Peak RAM vs. Number of Genes\n(sample count = %d; mean ± 95%% CI over %d runs)", n_samples_fixed, n_reps),
    x = "Number of genes",
    y = "Peak RAM (MiB)",
    tag = "D"
  ) +
  theme_minimal(base_size = 12)

#write.csv(results_genes_df, "~/Downloads/runtime_altering_genes_try_2.csv")