scripts/time_vs_C_cluster.R

In rchan26/exp4FusionRCPP: Hierarchical and Sequential Monte Carlo Fusion for exp(-(x^4)/2)

library(exp4FusionRCPP)

seed <- 408

set.seed(seed)
denominator <- 2:16
input_samples <- list()
standard_results <- list()
hier_results <- list()
seq_results <- list()

for (i in 1:length(denominator)) {
  print(denominator[i])
  input_samples[[i]] <- base_rejection_sampler_exp_4(beta = 1/denominator[i], nsamples = 50000, proposal_mean = 0,
                                                     proposal_sd = 1.5, dominating_M = 1.75)

  curve(tempered_target_density_exp_4(x, beta = 1/denominator[i]), -4, 4,
        main = denominator[i], ylab = 'tempered pdf')
  for (j in 1:length(input_samples[[i]])) {
    lines(density(input_samples[[i]][[j]]), col = 'blue')
  }

  # standard fork and join
  print('performing standard MC fusion')
  fnj_fused <- parallel_h_fusion_exp_4_rcpp(N_schedule = 10000, mean = 0, time_schedule = 0.5,
                                            start_beta = 1/denominator[i], C_schedule = denominator[i],
                                            L = 2, base_samples = input_samples[[i]], seed = seed)

  standard_results[[i]] <- list(fnj_fused$time, fnj_fused$rho_acc, fnj_fused$Q_acc, fnj_fused$rhoQ_acc)

  # hierarchical if denominator[i] is 2, 4, 8, or 16
  if (denominator[i]==2) {
    print('performing hierarchical MC fusion')
    hier_fused <- parallel_h_fusion_exp_4_rcpp(N_schedule = 10000, mean = 0, time_schedule = 0.5,
                                               start_beta = 1/2, C_schedule = 2,
                                               L = 2, base_samples = input_samples[[i]],
                                               seed = seed)
  } else if (denominator[i]==4) {
    print('performing hierarchical MC fusion')
    hier_fused <- parallel_h_fusion_exp_4_rcpp(N_schedule = rep(10000, 2), mean = 0, time_schedule = rep(0.5, 2),
                                               start_beta = 1/4, C_schedule = rep(2, 2),
                                               L = 3, base_samples = input_samples[[i]],
                                               seed = seed)
  } else if (denominator[i]==8) {
    print('performing hierarchical MC fusion')
    hier_fused <- parallel_h_fusion_exp_4_rcpp(N_schedule = rep(10000, 3), mean = 0, time_schedule = rep(0.5, 3),
                                               start_beta = 1/8, C_schedule = rep(2, 3),
                                               L = 4, base_samples = input_samples[[i]],
                                               seed = seed)
  } else if (denominator[i]==16) {
    print('performing hierarchical MC fusion')
    hier_fused <- parallel_h_fusion_exp_4_rcpp(N_schedule = rep(10000, 4), mean = 0, time_schedule = rep(0.5, 4),
                                               start_beta = 1/16, C_schedule = rep(2, 4),
                                               L = 5, base_samples = input_samples[[i]],
                                               seed = seed)
  }

  if (denominator[i] %in% c(2, 4, 8, 16)) {
    hier_results[[i]] <- list('time' = hier_fused$time,
                              'overall_rho' = hier_fused$overall_rho,
                              'overall_Q' = hier_fused$overall_Q,
                              'overall_rhoQ' = hier_fused$overall_rhoQ,
                              'rho_per_level' = hier_fused$rho_acc,
                              'Q_per_level' = hier_fused$Q_acc,
                              'rhoQ_per_level' = hier_fused$rhoQ_acc)
  } else {
    hier_results[[i]] <- NULL
  }


  # sequential
  print('performing sequential MC fusion')
  seq_fused <- parallel_sequential_fusion_exp_4_rcpp(N_schedule = rep(10000, denominator[i]-1), mean = 0,
                                                     time_schedule = rep(0.5, denominator[i]-1),
                                                     start_beta = 1/denominator[i], base_samples = input_samples[[i]],
                                                     seed = seed)

  seq_results[[i]] <- list('time' = seq_fused$time,
                           'overall_rho' = seq_fused$overall_rho,
                           'overall_Q' = seq_fused$overall_Q,
                           'overall_rhoQ' = seq_fused$overall_rhoQ,
                           'rho_per_level' = seq_fused$rho_acc,
                           'Q_per_level' = seq_fused$Q_acc,
                           'rhoQ_per_level' = seq_fused$rhoQ_acc)

  curve(target_density_exp_4(x, mean = 0), -4, 4, ylim = c(0, 0.5), main = denominator[i])
  # lines(density(fnj_fused$samples[[1]]), col = 'orange')
  if (!is.null(hier_results[[i]])) {
    lines(density(hier_fused$samples[[1]]), col = 'blue')
  }
  lines(density(seq_fused$samples[[1]]), col = 'green')
}

######################################## running time

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[1]]), ylim = c(0, 40000),
     ylab = 'Running time in seconds', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[1]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(sum(hier_results[[1]]$time), sum(hier_results[[3]]$time),
                                 sum(hier_results[[7]]$time), sum(hier_results[[15]]$time)), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(sum(hier_results[[1]]$time), sum(hier_results[[3]]$time),
                                sum(hier_results[[7]]$time), sum(hier_results[[15]]$time)), col = 'blue')
points(x = 2:16, y = sapply(1:15, function(i) sum(seq_results[[i]]$time)), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) sum(seq_results[[i]]$time)), col = 'green')

#################### log

plot(x = 2:16, y = sapply(1:15, function(i) log(standard_results[[i]][[1]])), ylim = c(-1, 12),
     ylab = '(logarithm) Running time in seconds', xlab = 'Number of Subposteriors (C)', col = 'red', pch = 4, lwd = 1.5)
lines(x = 2:16, y = sapply(1:15, function(i) log(standard_results[[i]][[1]])), col = 'red')
points(x = c(2, 4, 8, 16), y = log(c(sum(hier_results[[1]]$time), sum(hier_results[[3]]$time),
                                     sum(hier_results[[7]]$time), sum(hier_results[[15]]$time))), col = 'blue', pch = 0, lwd = 1.5)
lines(x = c(2, 4, 8, 16), y = log(c(sum(hier_results[[1]]$time), sum(hier_results[[3]]$time),
                                    sum(hier_results[[7]]$time), sum(hier_results[[15]]$time))), col = 'blue', lty = 2)
points(x = 2:16, y = sapply(1:15, function(i) log(sum(seq_results[[i]]$time))), col = 'darkgreen', pch = 1, lwd = 1.5)
lines(x = 2:16, y = sapply(1:15, function(i) log(sum(seq_results[[i]]$time))), col = 'darkgreen', lty = 3)
legend(x = 2, y = 12, legend = c('standard', 'hierarchical', 'sequential'),
       lty = c(1, 2, 3), pch = c(4, 0, 1), col = c('red', 'blue', 'darkgreen'))

######################################## rho acceptance (overall)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[2]]), ylim = c(0, 1),
     ylab = 'Rho Acceptance Rate', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[2]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_rho, hier_results[[3]]$overall_rho,
                                 hier_results[[7]]$overall_rho, hier_results[[15]]$overall_rho), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_rho, hier_results[[3]]$overall_rho,
                                hier_results[[7]]$overall_rho, hier_results[[15]]$overall_rho), col = 'blue')
points(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_rho), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_rho), col = 'green')

######################################## Q acceptance (overall)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[3]]), ylim = c(0, 1),
     ylab = 'Q Acceptance Rate', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[3]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_Q, hier_results[[3]]$overall_Q,
                                 hier_results[[7]]$overall_Q, hier_results[[15]]$overall_Q), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_Q, hier_results[[3]]$overall_Q,
                                hier_results[[7]]$overall_Q, hier_results[[15]]$overall_Q), col = 'blue')
points(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_Q), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_Q), col = 'green')

######################################## rho.Q acceptance (overall)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[4]]), ylim = c(0, 1),
     ylab = 'Overall Acceptance Rate', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[4]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_rhoQ, hier_results[[3]]$overall_rhoQ,
                                 hier_results[[7]]$overall_rhoQ, hier_results[[15]]$overall_rhoQ), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(hier_results[[1]]$overall_rhoQ, hier_results[[3]]$overall_rhoQ,
                                hier_results[[7]]$overall_rhoQ, hier_results[[15]]$overall_rhoQ), col = 'blue')
points(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_rhoQ), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) seq_results[[i]]$overall_rhoQ), col = 'green')


########################################################################################################################

######################################## rho acceptance (mean)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[2]]), ylim = c(0, 1),
     ylab = 'Rho Acceptance Rate (mean)', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[2]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$rho_per_level), mean(hier_results[[3]]$rho_per_level),
                                 mean(hier_results[[7]]$rho_per_level), mean(hier_results[[15]]$rho_per_level)), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$rho_per_level), mean(hier_results[[3]]$rho_per_level),
                                mean(hier_results[[7]]$rho_per_level), mean(hier_results[[15]]$rho_per_level)), col = 'blue', lty = 2)
points(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$rho_per_level)), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$rho_per_level)), col = 'green', lty = 2)

######################################## Q acceptance (mean)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[3]]), ylim = c(0, 1),
     ylab = 'Q Acceptance Rate', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[3]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$Q_per_level), mean(hier_results[[3]]$Q_per_level),
                                 mean(hier_results[[7]]$Q_per_level), mean(hier_results[[15]]$Q_per_level)), col = 'blue', pch = 4)
lines(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$Q_per_level), mean(hier_results[[3]]$Q_per_level),
                                mean(hier_results[[7]]$Q_per_level), mean(hier_results[[15]]$Q_per_level)), col = 'blue', lty = 2)
points(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$Q_per_level)), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$Q_per_level)), col = 'green', lty = 2)

######################################## overall acceptance (mean)

plot(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[4]]), ylim = c(0, 1),
     ylab = 'Overall Acceptance Rate', xlab = 'Number of Subposteriors (C)', col = 'orange', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) standard_results[[i]][[4]]), col = 'orange')
points(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$rhoQ_per_level), mean(hier_results[[3]]$rhoQ_per_level),
                                 mean(hier_results[[7]]$rhoQ_per_level), mean(hier_results[[15]]$rhoQ_per_level)), col = 'blue',pch = 4)
lines(x = c(2, 4, 8, 16), y = c(mean(hier_results[[1]]$rhoQ_per_level), mean(hier_results[[3]]$rhoQ_per_level),
                                mean(hier_results[[7]]$rhoQ_per_level), mean(hier_results[[15]]$rhoQ_per_level)), col = 'blue', lty = 2)
points(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$rhoQ_per_level)), col = 'green', pch = 4)
lines(x = 2:16, y = sapply(1:15, function(i) mean(seq_results[[i]]$rhoQ_per_level)), col = 'green', lty = 2)