library(microbenchmark) library(tensorPredictors) setwd("~/Work/tensorPredictors/sim/") base.name <- "sim_ising_perft" # Number of replications, sufficient for the performance test reps <- 5 # Sets the dimensions to be tested for runtime per method configs <- list( exact = list( # Exact method dim = 1:24, use_MC = FALSE, nr_threads = 1L # ignored in this case, but no special case neded ), MC = list( # Monte-Carlo Estimate dim = c(1:20, (3:13) * 10), use_MC = TRUE, nr_threads = 1L # default nr. of threads ), MC8 = list( # Monte-Carlo Estimate using 8 threads dim = c(1:20, (3:13) * 10), use_MC = TRUE, nr_threads = 8L # my machines nr of (virtual) cores ) ) # Simple function creating a parameter vector to be passed to `ising_m2`, the values # are irrelevant while its own execution time is (more or less) neglectable params <- function(dim) double(dim * (dim + 1L) / 2L) # Build expressions to be past to `microbenchmark` for performance testing expressions <- Reduce(c, Map(function(method) { config <- configs[[method]] Map(function(dim) { as.call(list(quote(ising_m2), params = substitute(params(dim), list(dim = dim)), use_MC = config$use_MC, nr_threads = config$nr_threads)) }, config$dim) }, names(configs))) # Performance tests perft.results <- microbenchmark(list = expressions, times = reps) # Convert performance test results to data frame for further processing (perft <- summary(perft.results)) # Ploting the performance simulation ################################################################################ ### TODO: Fix plotting ### ################################################################################ if (FALSE) { with(sim, { par(mfrow = c(2, 2), mar = c(5, 4, 4, 4) + 0.1) # Effect of Nr. of samples plot(range(nr_samples), range(mse - sd, mse + sd), type = "n", bty = "n", log = "xy", yaxt = "n", xlab = "Nr. Samples", ylab = "MSE", main = "Sample Size Effect (MSE)") groups <- split(sim, warmup) for (i in seq_along(groups)) { with(groups[[i]], { lines(nr_samples, mse, col = i, lwd = 2, type = "b", pch = 16) lines(nr_samples, mse - sd, col = i, lty = 2) lines(nr_samples, mse + sd, col = i, lty = 2) }) } right <- nr_samples == max(nr_samples) axis(4, at = mse[right], labels = warmup[right], lwd = 0, lwd.ticks = 1, line = -1.5) mtext("Warmup", side = 4, line = 1.5, at = exp(mean(range(log(mse[right]))))) y.pow <- -10:-1 axis(2, at = c(1, 10^y.pow), labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # Effect warmup length plot(range(warmup + 1), range(mse - sd, mse + sd), type = "n", bty = "n", log = "xy", xaxt = "n", yaxt = "n", xlab = "Warmup", ylab = "MSE", main = "Warmup Effect (MSE)") axis(1, warmup + 1, labels = as.integer(warmup)) groups <- split(sim, nr_samples) for (i in seq_along(groups)) { with(groups[[i]], { lines(warmup + 1, mse, col = i, lwd = 2, type = "b", pch = 16) lines(warmup + 1, mse - sd, col = i, lty = 2) lines(warmup + 1, mse + sd, col = i, lty = 2) }) } right <- warmup == max(warmup) axis(4, at = mse[right], labels = nr_samples[right], lwd = 0, lwd.ticks = 1, line = -1.5) mtext("Nr. Samples", side = 4, line = 1.5, at = exp(mean(range(log(mse[right]))))) axis(2, at = c(1, 10^y.pow), labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # Effect of Nr. of samples plot(range(nr_samples), range(merr), type = "n", bty = "n", log = "xy", yaxt = "n", xlab = "Nr. Samples", ylab = "Max Abs Error Mean", main = "Sample Size Effect (Abs Error)") groups <- split(sim, warmup) for (i in seq_along(groups)) { with(groups[[i]], { lines(nr_samples, merr, col = i, lwd = 2, type = "b", pch = 16) }) } right <- nr_samples == max(nr_samples) axis(4, at = merr[right], labels = warmup[right], lwd = 0, lwd.ticks = 1, line = -1.5) mtext("Warmup", side = 4, line = 1.5, at = exp(mean(range(log(merr[right]))))) y.pow <- -10:-1 axis(2, at = c(1, 10^y.pow), labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # Effect of warmup length plot(range(warmup + 1), range(merr), type = "n", bty = "n", log = "xy", xaxt = "n", yaxt = "n", xlab = "Warmup", ylab = "Max Abs Error Mean", main = "Warmup Effect (Abs Error)") axis(1, warmup + 1, labels = as.integer(warmup)) groups <- split(sim, nr_samples) for (i in seq_along(groups)) { with(groups[[i]], { lines(warmup + 1, merr, col = i, lwd = 2, type = "b", pch = 16) }) } right <- warmup == max(warmup) axis(4, at = merr[right], labels = nr_samples[right], lwd = 0, lwd.ticks = 1, line = -1.5) mtext("Nr. Samples", side = 4, line = 1.5, at = exp(mean(range(log(merr[right]))))) axis(2, at = c(1, 10^y.pow), labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) }) # Add common title mtext(main, side = 3, line = -2, outer = TRUE, font = 2, cex = 1.5) } # test_unscaled_prob <- function() { # test <- function(p) { # y <- sample.int(2, p, replace = TRUE) - 1L # theta <- vech(matrix(rnorm(p^2), p)) # C <- ising_m2(y, theta) # R <- exp(sum(vech(outer(y, y)) * theta)) # if (all.equal(C, R) == TRUE) { # cat("\033[92mSUCCESS: ") # } else { # cat("\033[91mFAILED: ") # } # cat(sprintf("p = %d, C = %e, R = %e\n", p, C, R)) # cat(" ", paste0(format(seq_along(y) - 1), collapse = " "), "\n") # cat(" y = ", paste0(c(".", "1")[y + 1], collapse = " ")) # cat("\033[0m\n\n") # } # devtools::load_all() # for (p in c(1, 10, 30:35, 62:66, 70, 128, 130)) { # test(p) # } # } # test_ising_sample <- function() { # test <- function(p) { # # theta <- vech(matrix(rnorm(p^2), p)) # # theta <- vech(matrix(0, p, p)) # theta <- -0.01 * vech(1 - diag(p)) # # theta <- vech(0.2 * diag(p)) # sample <- ising_sample(11, theta) # print.table(sample, zero.print = ".") # print(mean(sample)) # } # devtools::load_all() # for (p in c(1, 10, 30:35, 62:66, 70, 128, 130)) { # test(p) # } # } # test_ising_partition_func_MC <- function() { # test <- function(p) { # # theta <- vech(matrix(rnorm(p^2), p)) # # theta <- vech(matrix(0, p, p)) # theta <- -0.01 * vech(1 - diag(p)) # time_gmlm <- system.time(val_gmlm <- ising_partition_func_MC(theta)) # time_gmlm <- round(1000 * time_gmlm[["elapsed"]]) # if (p < 21) { # time_mvb <- system.time(val_mvb <- 1 / mvbernoulli::ising_prob0(theta)) # time_mvb <- round(1000 * time_mvb[["elapsed"]]) # } else { # val_mvb <- NaN # time_mvb <- -1 # } # cat(sprintf( # "dim = %d\n GMLM: time = %4d ms, val = %.4e\n MVB: time = %4d ms, val = %.4e\n", # p, time_gmlm, val_gmlm, time_mvb, val_mvb)) # } # devtools::load_all() # system.time( # # for (p in c(1, 10, 20, 30:35, 64, 70, 128, 130)) { # for (p in c(1, 10, 20, 30:35, 64)) { # test(p) # } # ) # } # # test_ising_partition_func_MC() # validate_ising_partition_func_MC <- function(theta_func) { # est_var <- function(dim) { # theta <- theta_func(dim) # time <- system.time(rep <- replicate(100, ising_partition_func_MC(theta))) # cat(sprintf("dim = %d, time = %.2e s, mean = %.2e, std.dev = %.2e\n", # dim, time[["elapsed"]], mean(rep), sd(rep))) # } # for (dim in 10 * (1:13)) { # est_var(dim) # } # } # # validate_ising_partition_func_MC(function(dim) { vech(matrix(rnorm(dim^2), dim)) }) # # validate_ising_partition_func_MC(function(dim) { vech(matrix(0, dim, dim)) }) # # validate_ising_partition_func_MC(function(dim) { -0.01 * vech(1 - diag(dim)) }) # # validate_ising_partition_func_MC(function(dim) { vech(0.2 * diag(dim)) }) # test_ising_partition_func_exact <- function(theta_func) { # test <- function(dim) { # theta <- theta_func(dim) # reps <- if (dim < 10) 100 else 10 # time <- system.time(replicate(reps, ising_partition_func_exact(theta))) # time <- time[["elapsed"]] / reps # cat(sprintf("dim = %d, time = %.2e s\n", dim, time)) # } # for (dim in 1:20) { # test(dim) # } # } # test_ising_partition_func_exact(function(dim) { vech(matrix(rnorm(dim^2), dim)) }) # ### Performance Measurement/Comparison # local({ # perft_exact <- local({ # dims <- 2:22 # cat("Exact perft:\n") # times <- sapply(dims, function(dim) { # reps <- if (dim < 10) 1000 else if (dim < 15) 100 else if (dim < 20) 10 else 4 # theta <- vech(matrix(rnorm(dim^2), dim)) # time <- system.time(replicate(reps, # ising_m2(theta, use_MC = FALSE) # )) # time <- time[["elapsed"]] / reps # cat(sprintf(" dim = %3d, reps = %3d, time per rep = %.2e s\n", dim, reps, time)) # time # }) # list(dims = dims, times = times) # }) # perft_MC <- local({ # dims <- c(2:21, 30, 40, 70, 100) # cat("Monte-Carlo perft:\n") # times <- sapply(dims, function(dim) { # reps <- if (dim < 20) 25 else if (dim < 40) 10 else 4 # theta <- vech(matrix(rnorm(dim^2), dim)) # time <- system.time(replicate(reps, # ising_m2(theta, use_MC = TRUE) # )) # time <- time[["elapsed"]] / reps # cat(sprintf(" dim = %3d, reps = %3d, time per rep = %.2e s\n", dim, reps, time)) # time # }) # list(dims = dims, times = times) # }) # perft_MC_thrd <- local({ # dims <- c(2:21, 30, 40, 70, 100) # cat("Monte-Carlo Multi-Threaded perft:\n") # times <- sapply(dims, function(dim) { # reps <- if (dim < 15) 25 else if (dim < 40) 10 else 4 # theta <- vech(matrix(rnorm(dim^2), dim)) # time <- system.time(replicate(reps, # ising_m2(theta, use_MC = TRUE, nr_threads = 6L) # )) # time <- time[["elapsed"]] / reps # cat(sprintf(" dim = %3d, reps = %3d, time per rep = %.2e s\n", dim, reps, time)) # time # }) # list(dims = dims, times = times) # }) # # Plot results # par(mfrow = c(1, 1)) # plot( # range(c(perft_MC_thrd$dims, perft_MC$dims, perft_exact$dims)), # range(c(perft_MC_thrd$times, perft_MC$times, perft_exact$times)), # type = "n", log = "xy", xlab = "Dimension p", ylab = "Time", xaxt = "n", yaxt = "n", # main = "Performance Comparison" # ) # # Add custom Y-axis # x.major.ticks <- as.vector(outer(c(2, 5, 10), 10^(0:5))) # x.minor.ticks <- as.vector(outer(2:10, 10^(0:5))) # axis(1, x.major.ticks, labels = as.integer(x.major.ticks)) # axis(1, x.minor.ticks, labels = NA, tcl = -0.25, lwd = 0, lwd.ticks = 1) # abline(v = x.major.ticks, col = "gray", lty = "dashed") # abline(v = x.minor.ticks, col = "lightgray", lty = "dotted") # # Add custom Y-axis # y.major.ticks <- c(10^(-9:1), 60, 600, 3600) # y.labels <- c( # expression(paste(n, s)), # expression(paste(10, n, s)), # expression(paste(100, n, s)), # expression(paste(mu, s)), # expression(paste(10, mu, s)), # expression(paste(100, mu, s)), # expression(paste(1, m, s)), # expression(paste(10, m, s)), # expression(paste(100, m, s)), # expression(paste(1, s)), # expression(paste(10, s)), # expression(paste(1, min)), # expression(paste(10, min)), # expression(paste(1, h)) # ) # y.minor.ticks <- c(as.vector(outer((1:10), 10^(-10:0))), 10 * (1:6), 60 * (2:10), 600 * (2:6)) # axis(2, at = y.major.ticks, labels = y.labels) # axis(2, at = y.minor.ticks, labels = NA, tcl = -0.25, lwd = 0, lwd.ticks = 1) # abline(h = y.major.ticks, col = "gray", lty = "dashed") # abline(h = y.minor.ticks, col = "lightgray", lty = "dotted") # legend("bottomright", col = c("red", "darkgreen", "blue"), lty = c(1, 1, 1), # bg = "white", # legend = c( # expression(paste("Exact ", O(2^p))), # expression(paste("MC ", O(p^2))), # expression(paste("MC Thrd ", O(p^2))) # ) # ) # with(perft_exact, { # points(dims, times, pch = 16, col = "red") # with(list(dims = tail(dims, -4), times = tail(times, -4)), { # lines(dims, exp(predict(lm(log(times) ~ dims))), col = "red") # }) # }) # with(perft_MC, { # points(dims, times, pch = 16, col = "darkgreen") # lines(dims, predict(lm(sqrt(times) ~ dims))^2, col = "darkgreen") # }) # with(perft_MC_thrd, { # points(dims, times, pch = 16, col = "blue") # lines(dims, predict(lm(sqrt(times) ~ dims))^2, col = "blue") # }) # }) # # dim <- 10 # # theta <- vech(matrix(rnorm(dim^2, 0, 1), dim, dim)) # # nr_threads <- 6L # # (m2.exact <- ising_m2(theta, use_MC = FALSE)) # # (m2.MC <- ising_m2(theta, use_MC = TRUE)) # # (m2.MC_thrd <- ising_m2(theta, use_MC = TRUE, nr_threads = nr_threads)) # # tcrossprod(ising_sample(1e4, theta)) / 1e4 # local({ # dim <- 20 # theta <- vech(matrix(rnorm(dim^2, 0, 1), dim, dim)) # A <- matrix(NA_real_, dim, dim) # A[lower.tri(A, diag = TRUE)] <- theta # A[lower.tri(A)] <- A[lower.tri(A)] / 2 # A[upper.tri(A)] <- t(A)[upper.tri(A)] # nr_threads <- 6L # time.exact <- system.time(m2.exact <- ising_m2(theta, use_MC = FALSE)) # time.MC <- system.time(m2.MC <- ising_m2(theta, use_MC = TRUE)) # time.MC_thrd <- system.time(m2.MC_thrd <- ising_m2(A, use_MC = TRUE, nr_threads = nr_threads)) # time.sample <- system.time(m2.sample <- tcrossprod(ising_sample(1e4, theta)) / 1e4) # range <- range(m2.exact, m2.MC, m2.MC_thrd) # par(mfrow = c(2, 2)) # matrixImage(m2.exact, main = sprintf("M2 Exact (time %.2f s)", time.exact[["elapsed"]]), zlim = range) # matrixImage(m2.MC, main = sprintf("M2 MC (time %.2f s)", time.MC[["elapsed"]]), zlim = range) # matrixImage(m2.MC_thrd, main = sprintf("M2 MC (%d threads, time %.2f s)", nr_threads, time.MC_thrd[["elapsed"]]), zlim = range) # matrixImage(m2.sample, main = sprintf("E_n(X X') (time %.2f s)", time.sample[["elapsed"]]), zlim = range) # # matrixImage(abs(m2.exact - m2.MC), main = "Abs. Error (Exact to MC)", zlim = c(-1, 1)) # }) # # Simulation # dims <- c(5, 10, 15, 20) # config.grid <- expand.grid( # nr_samples = c(10, 100, 1000, 10000), # warmup = c(0, 2, 10, 100), # dim = dims # ) # params <- Map(function(dim) vech(matrix(rnorm(dim^2, 0, 1), dim, dim)), dims) # names(params) <- dims # m2s.exact <- Map(ising_m2, params, use_MC = FALSE) # sim <- data.frame(t(apply(config.grid, 1, function(conf) { # # get same theta for every dimension # theta <- params[[as.character(conf["dim"])]] # m2.exact <- m2s.exact[[as.character(conf["dim"])]] # rep <- replicate(25, { # time <- system.time( # m2.MC <- ising_m2(theta, nr_samples = conf["nr_samples"], warmup = conf["warmup"], use_MC = TRUE) # ) # c(mse = mean((m2.exact - m2.MC)^2), err = max(abs(m2.exact - m2.MC)), time = time[["elapsed"]]) # }) # cat(sprintf("dim = %d, nr_samples = %6d, warmup = %6d, mse = %.4f\n", # conf["dim"], conf["nr_samples"], conf["warmup"], mean(rep["mse", ]))) # c( # conf, # mse = mean(rep["mse", ]), sd = sd(rep["mse", ]), merr = mean(rep["err", ]), # time = mean(rep["time", ]) # ) # }))) # plot.sim <- function(sim, main) { # with(sim, { # par(mfrow = c(2, 2), mar = c(5, 4, 4, 4) + 0.1) # # Effect of Nr. of samples # plot(range(nr_samples), range(mse - sd, mse + sd), # type = "n", bty = "n", log = "xy", yaxt = "n", # xlab = "Nr. Samples", ylab = "MSE", # main = "Sample Size Effect (MSE)") # groups <- split(sim, warmup) # for (i in seq_along(groups)) { # with(groups[[i]], { # lines(nr_samples, mse, col = i, lwd = 2, type = "b", pch = 16) # lines(nr_samples, mse - sd, col = i, lty = 2) # lines(nr_samples, mse + sd, col = i, lty = 2) # }) # } # right <- nr_samples == max(nr_samples) # axis(4, at = mse[right], labels = warmup[right], lwd = 0, lwd.ticks = 1, line = -1.5) # mtext("Warmup", side = 4, line = 1.5, at = exp(mean(range(log(mse[right]))))) # y.pow <- -10:-1 # axis(2, at = c(1, 10^y.pow), # labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # # Effect warmup length # plot(range(warmup + 1), range(mse - sd, mse + sd), # type = "n", bty = "n", log = "xy", xaxt = "n", yaxt = "n", # xlab = "Warmup", ylab = "MSE", # main = "Warmup Effect (MSE)") # axis(1, warmup + 1, labels = as.integer(warmup)) # groups <- split(sim, nr_samples) # for (i in seq_along(groups)) { # with(groups[[i]], { # lines(warmup + 1, mse, col = i, lwd = 2, type = "b", pch = 16) # lines(warmup + 1, mse - sd, col = i, lty = 2) # lines(warmup + 1, mse + sd, col = i, lty = 2) # }) # } # right <- warmup == max(warmup) # axis(4, at = mse[right], labels = nr_samples[right], lwd = 0, lwd.ticks = 1, line = -1.5) # mtext("Nr. Samples", side = 4, line = 1.5, at = exp(mean(range(log(mse[right]))))) # axis(2, at = c(1, 10^y.pow), # labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # # Effect of Nr. of samples # plot(range(nr_samples), range(merr), # type = "n", bty = "n", log = "xy", yaxt = "n", # xlab = "Nr. Samples", ylab = "Max Abs Error Mean", # main = "Sample Size Effect (Abs Error)") # groups <- split(sim, warmup) # for (i in seq_along(groups)) { # with(groups[[i]], { # lines(nr_samples, merr, col = i, lwd = 2, type = "b", pch = 16) # }) # } # right <- nr_samples == max(nr_samples) # axis(4, at = merr[right], labels = warmup[right], lwd = 0, lwd.ticks = 1, line = -1.5) # mtext("Warmup", side = 4, line = 1.5, at = exp(mean(range(log(merr[right]))))) # y.pow <- -10:-1 # axis(2, at = c(1, 10^y.pow), # labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # # Effect of warmup length # plot(range(warmup + 1), range(merr), # type = "n", bty = "n", log = "xy", xaxt = "n", yaxt = "n", # xlab = "Warmup", ylab = "Max Abs Error Mean", # main = "Warmup Effect (Abs Error)") # axis(1, warmup + 1, labels = as.integer(warmup)) # groups <- split(sim, nr_samples) # for (i in seq_along(groups)) { # with(groups[[i]], { # lines(warmup + 1, merr, col = i, lwd = 2, type = "b", pch = 16) # }) # } # right <- warmup == max(warmup) # axis(4, at = merr[right], labels = nr_samples[right], lwd = 0, lwd.ticks = 1, line = -1.5) # mtext("Nr. Samples", side = 4, line = 1.5, at = exp(mean(range(log(merr[right]))))) # axis(2, at = c(1, 10^y.pow), # labels = c(1, sapply(y.pow, function(pow) eval(substitute(expression(10^i), list(i = pow)))))) # }) # # Add common title # mtext(main, side = 3, line = -2, outer = TRUE, font = 2, cex = 1.5) # } # plot.sim(subset(sim, sim$dim == 5), main = "Dim = 5") # plot.sim(subset(sim, sim$dim == 10), main = "Dim = 10") # plot.sim(subset(sim, sim$dim == 15), main = "Dim = 15") # plot.sim(subset(sim, sim$dim == 20), main = "Dim = 20") # dim <- 120 # params <- rnorm(dim * (dim + 1) / 2) # A <- matrix(NA_real_, dim, dim) # A[lower.tri(A, diag = TRUE)] <- params # A[lower.tri(A)] <- A[lower.tri(A)] / 2 # A[upper.tri(A)] <- t(A)[upper.tri(A)] # seed <- abs(as.integer(100000 * rnorm(1))) # all.equal( # { set.seed(seed); ising_sample(11, params) }, # { set.seed(seed); ising_sample(11, A) } # ) # x <- sample(0:1, 10, TRUE) # sum(vech(outer(x, x)) * params) # sum(x * (A %*% x)) # # M <- matrix(NA, dim, dim) # # M[lower.tri(M, diag = TRUE)] <- seq_len(dim * (dim + 1) / 2) - 1 # # rownames(M) <- (1:dim) - 1 # # colnames(M) <- (1:dim) - 1 # # print.table(M) # # i <- seq(0, dim - 1) # # (i * (2 * dim + 1 - i)) / 2 # # I <- 0 # # for (i in seq(0, dim - 1)) { # # print(I) # # I <- I + dim - i # # } # m2.exact <- vech.pinv(ising_m2(params, use_MC = FALSE)) # m2.MC <- vech.pinv(ising_m2(params, use_MC = TRUE)) # m2.mat <- tcrossprod(ising_sample(1e4, A)) / 1e4 # m2.vech <- tcrossprod(ising_sample(1e4, params)) / 1e4 # par(mfrow = c(2, 2)) # matrixImage(m2.exact, main = "exact") # matrixImage(m2.MC, main = "MC") # matrixImage(m2.mat, main = "Sample mat") # matrixImage(m2.vech, main = "Sample vech")