scripts/hline_tests.R

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

parallel_line_fusion_exp_4 <- function(N_schedule, time_schedule, start_beta, C_schedule, L, base_samples, study = F) {
  if (length(N_schedule) != (L-1)) stop('length of N_schedule must be equal to (L-1)')
  if (length(time_schedule) != (L-1)) stop('length of time_schedule must be equal to (L-1)')
  if (length(C_schedule) == (L-1)) {
~    # check that at each level, we are fusing a suitable number
    for (l in length(C_schedule):1) {
      if (((1/start_beta)/prod(C_schedule[(L-1):l]))%%1 != 0) {
        stop('check that (1/beta)/prod(C_schedule[(L-1):l]) is an integer for l=L-1,...,1')
      }
    }
  } else {
    stop('C_schedule must be a vector of length (L-1)')
  }
  
  # we append 1 to the vector C_schedule to make the indices work later on when we call fusion
  # we need this so that we can set the right value for beta when fusing up the levels
  C_schedule <- c(C_schedule, 1)
  
  # initialising study results
  hier_samples <- list()
  hier_samples[[L]] <- base_samples # base level
  time_taken_per_level <- rep(0, L-1)
  input_betas <- list()
  input_betas[[L]] <- NA
  output_beta <- c(1:(L-1), start_beta)
  
  # make some vectors for acceptance rates of the level
  rho_acc <- rep(0, L-1)
  Q_acc <- rep(0, L-1)
  rhoQ_acc <- rep(0, L-1)
  
  # creating parallel cluster
  n_cores <- parallel::detectCores()
  cl <- parallel::makeCluster(n_cores, outfile = 'parallel_line_fusion_exp_4.txt')
  
  # creating variable and functions list to pass into cluster using clusterExport
  varlist <- list("phi_function_exp_4", "bounds_phi_function_exp_4", "simulate_langevin_diffusion_exp_4",
                  "fusion_diff_exp_4", "N_schedule", "time_schedule", "start_beta", "C_schedule", "L", "base_samples")
  parallel::clusterExport(cl, envir = environment(), varlist = varlist)
  
  # exporting functions from layeredBB package to simulate layered Brownian bridges
  parallel::clusterExport(cl, varlist = ls("package:layeredBB"))
  
  # add to output file that starting hierarchical fusion
  cat('Starting hierarchical fusion \n', file = 'parallel_line_fusion_exp_4.txt')
  
  # parallelising tasks for each level going up the hiearchy
  for (k in ((L-1):1)) {
    samples_per_core <- rep(floor(N_schedule[k]/n_cores), n_cores)
    if (sum(samples_per_core) != N_schedule[k]) {
      remainder <- N_schedule[k] %% n_cores
      samples_per_core[1:remainder] <- samples_per_core[1:remainder] + 1
    }
    
    # performing Fusion for this level
    # printing out some stuff to log file to track the progress
    cat('########################\n', file = 'parallel_line_fusion_exp_4.txt', append = T)
    cat('Starting to fuse', C_schedule[k], 'densities of pi^beta, where beta =', prod(C_schedule[L:(k+1)]), '/',
        (1/start_beta), 'for level', k, 'with time', time_schedule[k], ', which is using', n_cores, 'cores\n',
        file = 'parallel_line_fusion_exp_4.txt', append = T)
    cat('Obtaining', N_schedule[k], 'samples for beta =', prod(C_schedule[L:k]), '/', (1/start_beta),
        '\n', file = 'parallel_line_fusion_exp_4.txt', append = T)
    cat('########################\n', file = 'parallel_line_fusion_exp_4.txt', append = T)
    
    # starting fusion
    pcm <- proc.time()
    fused <- parallel::parLapplyLB(cl, X = 1:length(samples_per_core), fun = function(i) {
      fusion_diff_exp_4(N = samples_per_core[i],
                        time = time_schedule[k],
                        C = C_schedule[k],
                        samples_to_fuse = rep(list(hier_samples[[k+1]]), C_schedule[k]),
                        betas = prod(C_schedule[L:(k+1)])*(start_beta),
                        level = k,
                        acceptance_rate = TRUE)})
    final <- proc.time() - pcm
    
    # putting the fusion samples into hier_samples[[k]]
    hier_samples[[k]] <- unlist(sapply(1:length(fused), function(i) fused[[i]]$samples))
    
    if (study) {
      # calculating the acceptance rates for all nodes in the current level
      # summing all the iterations for each core
      rho_iterations <- sum(sapply(1:length(fused), function(i) fused[[i]]$iters_rho))
      Q_iterations <- sum(sapply(1:length(fused), function(i) fused[[i]]$iters_Q))
      
      rho_acc[k] <- Q_iterations / rho_iterations
      Q_acc[k] <- (N_schedule[k]) / Q_iterations
      rhoQ_acc[k] <- (N_schedule[k]) / rho_iterations
      time_taken_per_level[k] <- final['elapsed']
      input_betas[[k]] <- rep(prod(C_schedule[L:(k+1)])*(start_beta), C_schedule[k])
      output_beta[k] <- prod(C_schedule[L:k])*(start_beta)
    }
  }
  
  # stopping cluster
  parallel::stopCluster(cl)
  # print completion
  cat('Completed hierarchical fusion\n', file = 'parallel_line_fusion_exp_4.txt', append = T)
  
  if (study) {
    return(list('samples' = hier_samples, 'time' = time_taken_per_level,
                'rho_acc' = rho_acc, 'Q_acc' = Q_acc, 'rhoQ_acc' = rhoQ_acc,
                'input_betas' = input_betas, 'output_beta' = output_beta, 'diffusion_times' = time_schedule))
  } else {
    return(hier_samples)
  }
}

parallel_line_fusion_TA_exp_4 <- function(N_schedule, global_T, start_beta, C_schedule, L, base_samples, study = F) {
  if (length(N_schedule) != (L-1)) stop('length of N_schedule must be equal to (L-1)')
  if (length(C_schedule) == (L-1)) {
    # check that at each level, we are fusing a suitable number
    for (l in length(C_schedule):1) {
      if (((1/start_beta)/prod(C_schedule[(L-1):l]))%%1 != 0) {
        stop('check that (1/beta)/prod(C_schedule[(L-1):l]) is an integer for l=L-1,...,1')
      }
    }
  } else {
    stop('C_schedule must be a vector of length (L-1)')
  }
  
  # change global_T by multiplying it by sqrt(start_beta)
  global_T <- global_T*sqrt(start_beta)
  
  # we append 1 to the vector C_schedule to make the indices work later on when we call fusion
  # we need this so that we can set the right value for beta when fusing up the levels
  C_schedule <- c(C_schedule, 1)
  
  # initialising study results
  hier_samples <- list()
  hier_samples[[L]] <- base_samples # base level
  time_taken_per_level <- rep(0, L-1)
  input_betas <- list()
  input_betas[[L]] <- NA
  output_beta <- c(1:(L-1), start_beta)
  diffusion_times <- list()
  
  # make some vectors for acceptance rates of the level
  rho_acc <- rep(0, L-1)
  Q_acc <- rep(0, L-1)
  rhoQ_acc <- rep(0, L-1)
  
  # creating parallel cluster
  n_cores <- parallel::detectCores()
  cl <- parallel::makeCluster(n_cores, outfile = 'parallel_line_fusion_TA_exp_4.txt')
  
  # creating variable and functions list to pass into cluster using clusterExport
  varlist <- list("phi_function_exp_4", "bounds_phi_function_exp_4", "simulate_langevin_diffusion_exp_4",
                  "fusion_TA_exp_4", "N_schedule", "global_T", "start_beta", "C_schedule", "L", "base_samples")
  parallel::clusterExport(cl, envir = environment(), varlist = varlist)
  
  # exporting functions from layeredBB package to simulate layered Brownian bridges
  parallel::clusterExport(cl, varlist = ls("package:layeredBB"))
  
  # add to output file that starting hierarchical fusion
  cat('Starting hierarchical fusion \n', file = 'parallel_line_fusion_TA_exp_4.txt')
  
  # parallelising tasks for each level going up the hiearchy
  for (k in ((L-1):1)) {
    samples_per_core <- rep(floor(N_schedule[k]/n_cores), n_cores)
    if (sum(samples_per_core) != N_schedule[k]) {
      remainder <- N_schedule[k] %% n_cores
      samples_per_core[1:remainder] <- samples_per_core[1:remainder] + 1
    }
    
    # performing Fusion for this level
    # printing out some stuff to log file to track the progress
    cat('########################\n', file = 'parallel_line_fusion_TA_exp_4.txt', append = T)
    cat('Starting to fuse', C_schedule[k], 'densities of pi^beta, where beta =', prod(C_schedule[L:(k+1)]), '/', 
        (1/start_beta), 'for level', k, 'with time', global_T/sqrt(prod(C_schedule[L:(k+1)])*(start_beta)),
        ', which is using', n_cores, 'cores\n',
        file = 'parallel_line_fusion_TA_exp_4.txt', append = T)
    cat('Obtaining', N_schedule[k], 'samples for beta =', prod(C_schedule[L:k]), '/', (1/start_beta),
        '\n', file = 'parallel_line_fusion_TA_exp_4.txt', append = T)
    cat('########################\n', file = 'parallel_line_fusion_TA_exp_4.txt', append = T)
    
    # starting fusion
    pcm <- proc.time()
    fused <- parallel::parLapplyLB(cl, X = 1:length(samples_per_core), fun = function(i) {
      fusion_TA_exp_4(N = samples_per_core[i],
                      time = global_T,
                      C = C_schedule[k],
                      samples_to_fuse = rep(list(hier_samples[[k+1]]), C_schedule[k]),
                      betas = prod(C_schedule[L:(k+1)])*(start_beta),
                      sample_weights = sqrt(prod(C_schedule[L:(k+1)])*(start_beta)),
                      level = k,
                      acceptance_rate = TRUE)})
    final <- proc.time() - pcm
    
    # putting the fusion samples into hier_samples[[k]]
    hier_samples[[k]] <- unlist(sapply(1:length(fused), function(i) fused[[i]]$samples))
    
    if (study) {
      # calculating the acceptance rates for all nodes in the current level
      # summing all the iterations for each core
      rho_iterations <- sum(sapply(1:length(fused), function(i) fused[[i]]$iters_rho))
      Q_iterations <- sum(sapply(1:length(fused), function(i) fused[[i]]$iters_Q))
      
      rho_acc[k] <- Q_iterations / rho_iterations
      Q_acc[k] <- (N_schedule[k]) / Q_iterations
      rhoQ_acc[k] <- (N_schedule[k]) / rho_iterations
      time_taken_per_level[k] <- final['elapsed']
      input_betas[[k]] <- rep(prod(C_schedule[L:(k+1)])*(start_beta), C_schedule[k])
      output_beta[k] <- prod(C_schedule[L:k])*(start_beta)
      diffusion_times[[k]] <- global_T / rep(sqrt(prod(C_schedule[L:(k+1)])*(start_beta)), C_schedule[k])
    }
  }
  
  # stopping cluster
  parallel::stopCluster(cl)
  # print completion
  cat('Completed hierarchical fusion\n', file = 'parallel_line_fusion_TA_exp_4.txt', append = T)
  
  hier_samples[[1]] <- unlist(hier_samples[[1]])
  if (study) {
    return(list('samples' = hier_samples, 'time' = time_taken_per_level,
                'rho_acc' = rho_acc, 'Q_acc' = Q_acc, 'rhoQ_acc' = rhoQ_acc,
                'input_betas' = input_betas, 'output_beta' = output_beta, 'diffusion_times' = diffusion_times))
  } else {
    return(hier_samples)
  }
}


######################################## examples ########################################
library(exp4Tempering)

# obtaining samples for target via rejection sampling
test_target_mc <- sample_from_fc(N = 10000, proposal_mean = 0, proposal_sd = 1, dominating_M = 1.35, beta = 1)

######################################## beta = 1/4

# samples for base level (4 nodes for normal hierarchical fusion, just the first node (input_samples1[[1]]) for line fusion)
input_samples1 <- base_rejection_sampler_exp_4(beta = 1/4, nsamples = 100000,
                                               proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4)

# regular hierarchical fusion
test1_standard <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 2), time_schedule = rep(1, 2),
                                          start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
                                          base_samples = input_samples1, study = T)

test1_TA <- parallel_h_fusion_TA_exp_4(N_schedule = rep(10000, 2), global_T = 1,
                                       start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
                                       base_samples = input_samples1, study = T)

# line fusion
test1_line <- parallel_line_fusion_exp_4(N_schedule = rep(10000, 2), time_schedule = rep(1, 2),
                                          start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
                                          base_samples = input_samples1[[1]], study = T)

test1_line_TA <- parallel_line_fusion_TA_exp_4(N_schedule = rep(10000, 2), global_T = 1,
                                                 start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
                                                 base_samples = input_samples1[[1]], study = T)

# plot the densities from line fusion
curve(target_density_exp_4(x), -2.5, 2.5, ylim = c(0,1))
lines(density(test_target_mc), col = 'black')
lines(density(test1_line$samples[[1]]), col = 'green')
lines(density(test1_line_TA$samples[[1]]), col = 'blue')

# print differences in time between the standard approaches
print(test1_standard$time)
print(test1_line$time)

# print differences in time between the TA approaches
print(test1_TA$time)
print(test1_line_TA$time)

# print total variation from the samples to the true density
print(total_variation_true_exp_4(test_target_mc))
print(total_variation_true_exp_4(test1_standard$samples[[1]]))
print(total_variation_true_exp_4(test1_TA$samples[[1]]))
print(total_variation_true_exp_4(test1_line$samples[[1]]))
print(total_variation_true_exp_4(test1_line$samples[[1]]))

# acceptance rate plots
acceptance_rate_plots(hier1 = test1_standard, hier2 = test1_line, time = 1)
acceptance_rate_plots(hier1 = test1_TA, hier2 = test1_line_TA, time = 1)

######################################## beta = 1/8

# samples for base level (8 nodes for normal hierarchical fusion, just the first node (input_samples2[[1]]) for line fusion)
input_samples2 <- hmc_base_sampler_exp_4(beta = 1/8, nsamples = 100000, nchains = 8)

# regular hierarchical fusion
test2_standard <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 3), time_schedule = rep(1, 3),
                                          start_beta = 1/8, C_schedule = rep(2, 3), L = 4,
                                          base_samples = input_samples2, study = T)
test2_TA <- parallel_h_fusion_TA_exp_4(N_schedule = rep(10000, 3), global_T = 1,
                                       start_beta = 1/8, C_schedule = rep(2, 3), L = 4,
                                       base_samples = input_samples2, study = T)

# line fusion
test2_line <- parallel_line_fusion_exp_4(N_schedule = rep(10000, 3), time_schedule = rep(1, 3),
                                           start_beta = 1/8, C_schedule = rep(2, 3), L = 4,
                                           base_samples = input_samples2[[1]], study = T)

test2_line_TA <- parallel_line_fusion_TA_exp_4(N_schedule = rep(10000, 3), global_T = 1, 
                                                 start_beta = 1/8, C_schedule = rep(2, 3), L = 4, 
                                                 base_samples = input_samples2[[1]], study = T)

# print differences in time between the standard approaches
print(test2_standard$time)
print(test2_line$time)

# print differences in time between the TA approaches
print(test2_TA$time)
print(test2_line_TA$time)

# print total variation from the samples to the true density
print(total_variation_true_exp_4(test_target_mc))
print(total_variation_true_exp_4(test2_standard$samples[[1]]))
print(total_variation_true_exp_4(test2_TA$samples[[1]]))
print(total_variation_true_exp_4(test2_line$samples[[1]]))
print(total_variation_true_exp_4(test2_line$samples[[1]]))

# plot the densities from line fusion
curve(target_density_exp_4(x), -2.5, 2.5, ylim = c(0,1))
lines(density(test_target_mc), col = 'black')
lines(density(test2_line$samples[[1]]), col = 'green')
lines(density(test2_line_TA$samples[[1]]), col = 'blue')

# acceptance rate plots
acceptance_rate_plots(hier1 = test2_standard, hier2 = test2_line, time = 1)
acceptance_rate_plots(hier1 = test2_TA, hier2 = test2_line_TA, time = 1)