From 4324295162b622c84d066f1b357df8e16edcfcf8 Mon Sep 17 00:00:00 2001 From: daniel Date: Mon, 11 Dec 2023 13:23:29 +0100 Subject: [PATCH] changed isimg_m2() to use the LogSumExp trick to compute the log_prob_0 instead of prob_0, and som other stuff --- sim/sim_ising_perft.R | 610 ++++++++++++++++++++++ sim/sim_utils.R | 205 ++++++++ tensorPredictors/R/HOCCA.R | 55 ++ tensorPredictors/R/SIR.R | 36 ++ tensorPredictors/R/det.R | 7 + tensorPredictors/R/gmlm_ising.R | 19 +- tensorPredictors/R/ising_m2.R | 38 +- tensorPredictors/R/kronperm.R | 36 ++ tensorPredictors/R/mat_mani_projections.R | 120 +++++ tensorPredictors/R/pinv.R | 17 + tensorPredictors/R/slice_select.R | 67 +++ tensorPredictors/R/solve.R | 11 + tensorPredictors/R/sym.R | 81 +++ tensorPredictors/src/det.c | 73 +++ tensorPredictors/src/det.h | 15 + tensorPredictors/src/ising_m2.c | 68 ++- tensorPredictors/src/solve.c | 100 ++++ tensorPredictors/src/solve.h | 19 + 18 files changed, 1511 insertions(+), 66 deletions(-) create mode 100644 sim/sim_ising_perft.R create mode 100644 sim/sim_utils.R create mode 100644 tensorPredictors/R/HOCCA.R create mode 100644 tensorPredictors/R/SIR.R create mode 100644 tensorPredictors/R/det.R create mode 100644 tensorPredictors/R/kronperm.R create mode 100644 tensorPredictors/R/mat_mani_projections.R create mode 100644 tensorPredictors/R/pinv.R create mode 100644 tensorPredictors/R/slice_select.R create mode 100644 tensorPredictors/R/solve.R create mode 100644 tensorPredictors/R/sym.R create mode 100644 tensorPredictors/src/det.c create mode 100644 tensorPredictors/src/det.h create mode 100644 tensorPredictors/src/solve.c create mode 100644 tensorPredictors/src/solve.h diff --git a/sim/sim_ising_perft.R b/sim/sim_ising_perft.R new file mode 100644 index 0000000..ff7232b --- /dev/null +++ b/sim/sim_ising_perft.R @@ -0,0 +1,610 @@ +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") diff --git a/sim/sim_utils.R b/sim/sim_utils.R new file mode 100644 index 0000000..072d096 --- /dev/null +++ b/sim/sim_utils.R @@ -0,0 +1,205 @@ +#' Some utility function used in simulations + +#' Computes the orthogonal projection matrix onto the span of `A` +proj <- function(A) tcrossprod(A, A %*% solve(crossprod(A, A))) + + +#' Logging Errors and Warnings +#' +#' Register a global warning and error handler for logging warnings/errors with +#' current simulation repetition session informatin allowing to reproduce problems +#' +#' @examples +#' # Usage +#' globalCallingHandlers(list( +#' message = exceptionLogger("warning.log"), +#' warning = exceptionLogger("warning.log"), +#' error = exceptionLogger("error.log") +#' )) +#' # Do some stuff where an error might occure +#' for (rep in 1:1000) { +#' # Store additional information logged with an error when an exception occures +#' storeExceptionInfo(rep = rep) +#' # Do work +#' stopifnot(rep < 17) +#' } +#' +assign(".exception.info", NULL, env = .GlobalEnv) +exceptionLogger <- function(file.name) { + force(file.name) + function(ex) { + log <- file(file.name, open = "a+") + cat("\n### Log At: ", format(Sys.time()), "\n", file = log) + cat("# Exception:\n", file = log) + cat(as.character.error(ex), file = log) + cat("\n# Exception Info:\n", file = log) + dump(".exception.info", envir = .GlobalEnv, file = log) + cat("\n# Traceback:\n", file = log) + + # add Traceback (see: `traceback()` which the following is addapted from) + n <- length(x <- .traceback(NULL, max.lines = -1L)) + if (n == 0L) { + cat("No traceback available", "\n", file = log) + } else { + for (i in 1L:n) { + xi <- x[[i]] + label <- paste0(n - i + 1L, ": ") + m <- length(xi) + srcloc <- if (!is.null(srcref <- attr(xi, "srcref"))) { + srcfile <- attr(srcref, "srcfile") + paste0(" at ", basename(srcfile$filename), "#", srcref[1L]) + } + if (isTRUE(attr(xi, "truncated"))) { + xi <- c(xi, " ...") + m <- length(xi) + } + if (!is.null(srcloc)) { + xi[m] <- paste0(xi[m], srcloc) + } + if (m > 1) { + label <- c(label, rep(substr(" ", 1L, + nchar(label, type = "w")), m - 1L)) + } + cat(paste0(label, xi), sep = "\n", file = log) + } + } + close(log) + } +} + +#' Used in conjuntion with `exceptionLogger()` +storeExceptionInfo <- function(...) { + info <- list(...) + info$RNGking <- RNGkind() + info$.Random.seed <- get0(".Random.seed", envir = .GlobalEnv) + .GlobalEnv$.exception.info <- info +} + + +### Simulation logging routine +#' @examples +#' # Create a CSV logger +#' logger <- CSV.logger("test.csv", header = c("A", "B", "C")) +#' +#' # Store some values in current environment +#' A <- 1 +#' B <- 2 +#' # write values to csv file (order irrelevent) +#' logger(C = -3, B = -2) +#' +#' # another line +#' A <- 10 +#' B <- 20 +#' logger(B = -20, C = -30) +#' +#' # read the file back in +#' read.csv("test.csv") +#' +#' # In simulations it is often usefull to time stamp the files +#' nr <- 5 +#' logger <- CSV.logger( +#' sprintf("test-nr%03d-%s.csv", nr, format(Sys.time(), "%Y%m%dT%H%M")), +#' header = c("A", "B", "C") +#' ) +#' +CSV.logger <- function(file.name, header) { + force(file.name) + + # CSV header, used to ensure correct value/column mapping when writing to file + force(header) + cat(paste0(header, collapse = ","), "\n", sep = "", file = file.name) + + function(...) { + # get directly provided data + arg.data <- list(...) + # all arguments must be given with a name + if (length(arg.data) && is.null(names(arg.data))) { + stop("Arguments must be given with names") + } + # check if all elements have a described CSV header column + unknown <- !(names(arg.data) %in% header) + if (any(unknown)) { + stop("Got unknown columns: ", paste0(names(arg.data)[unknown], collapse = ", ")) + } + # get missing values from environment + missing <- !(header %in% names(arg.data)) + env <- parent.frame() + data <- c(arg.data, mget(header[missing], envir = env)) + # Format all aguments + data <- Map(format, data) + # collaps into single line + line <- paste0(data[header], collapse = ",") + # write data line to file + cat(line, "\n", sep = "", file = file.name, append = TRUE) + } +} + + +# Set colors for every method +methods <- c("gmlm", "pca", "hopca", "tsir", "mgcca", "lpca", "clpca", "tnormal", "sir") +col.methods <- palette.colors(n = length(methods), palette = "Okabe-Ito", recycle = FALSE) +names(col.methods) <- methods + + +# Comparison plot of one measure for a simulation +plot.sim <- function(sim, measure.name, ..., ylim = c(0, 1)) { + par.default <- par(pch = 16, bty = "n", lty = "solid", lwd = 1.5) + + # # Set colors for every method + # methods <- c("gmlm", "pca", "hopca", "tsir", "mgcca", "lpca", "clpca", "tnormal") + # col.methods <- palette.colors(n = length(methods), palette = "Okabe-Ito", recycle = FALSE) + # names(col.methods) <- methods + + # Remain sample size grouping variable to avoid conflicts + aggr.mean <- aggregate(sim, list(sampleSize = sim$sample.size), mean) + aggr.median <- aggregate(sim, list(sampleSize = sim$sample.size), median) + aggr.sd <- aggregate(sim, list(sampleSize = sim$sample.size), sd) + aggr.min <- aggregate(sim, list(sampleSize = sim$sample.size), min) + aggr.max <- aggregate(sim, list(sampleSize = sim$sample.size), max) + + with(aggr.mean, { + plot(range(sampleSize), ylim, type = "n", ...) + for (dist.name in ls(pattern = paste0("^", measure.name))) { + mean <- get(dist.name) + median <- aggr.median[aggr.sd$sampleSize == sampleSize, dist.name] + sd <- aggr.sd[aggr.sd$sampleSize == sampleSize, dist.name] + min <- aggr.min[aggr.sd$sampleSize == sampleSize, dist.name] + max <- aggr.max[aggr.sd$sampleSize == sampleSize, dist.name] + method <- tail(strsplit(dist.name, ".", fixed = TRUE)[[1]], 1) + col <- col.methods[method] + lines(sampleSize, mean, type = "o", col = col, lty = 1, lwd = 2 + (method == "gmlm")) + lines(sampleSize, mean + sd, col = col, lty = 2, lwd = 0.8) + lines(sampleSize, mean - sd, col = col, lty = 2, lwd = 0.8) + lines(sampleSize, median, col = col, lty = 1, lwd = 1) + lines(sampleSize, min, col = col, lty = 3, lwd = 0.6) + lines(sampleSize, max, col = col, lty = 3, lwd = 0.6) + } + + legend("topright", col = col.methods, lty = 1, legend = names(col.methods), + bty = "n", lwd = par("lwd"), pch = par("pch")) + }) + + # reset plotting default prameters + par(par.default) +} + + +timer.env <- new.env() +start.timer <- function() { + assign("start.time", proc.time()[["elapsed"]], envir = timer.env) +} +clear.timer <- function() { + assign("total.time", 0, envir = timer.env) +} +end.timer <- function() { + end.time <- proc.time()[["elapsed"]] + start.time <- get("start.time", envir = timer.env) + total.time <- get0("total.time", envir = timer.env) + if (is.null(total.time)) { + total.time <- 0 + } + elapsed <- end.time - start.time + total.time <- total.time + elapsed + assign("total.time", total.time, envir = timer.env) + c(elapsed = elapsed, total.time = total.time) +} diff --git a/tensorPredictors/R/HOCCA.R b/tensorPredictors/R/HOCCA.R new file mode 100644 index 0000000..c67efec --- /dev/null +++ b/tensorPredictors/R/HOCCA.R @@ -0,0 +1,55 @@ +#' A simple higher order (multi-way) canonical correlation analysis. +#' +#' @param X multi-dimensional array +#' @param Y multi-dimensional array with the same nr. of dimensions and equal +#' sample axis to `X`. +#' @param sample.axis integer indicationg which axis enumerates observations +#' +#' @export +HOCCA <- function(X, Y, sample.axis = length(dim(X)), centerX = TRUE, centerY = TRUE) { + + # ensure sample axis is the last axis + if (!missing(sample.axis)) { + modes <- seq_along(dim(X))[-sample.axis] + X <- aperm(X, c(modes, sample.axis)) + Y <- aperm(Y, c(modes, sample.axis)) + } + modes <- seq_len(length(dim(X)) - 1L) + dimX <- head(dim(X), -1L) + dimF <- head(dim(F), -1L) + sample.size <- tail(dim(X), 1L) + + # center `X` and `Y` + if (centerX) { + X <- X - as.vector(rowMeans(X, dims = length(dim(X)) - 1L)) + } + if (centerY) { + Y <- Y - as.vector(rowMeans(Y, dims = length(dim(Y)) - 1L)) + } + + # estimate marginal covariance matrices + CovXX <- Map(function(mode) mcrossprod(X, mode = mode) / prod(dim(X)[-mode]), modes) + CovYY <- Map(function(mode) mcrossprod(Y, mode = mode) / prod(dim(Y)[-mode]), modes) + # and the "covariance tensor" + CovXY <- array(tcrossprod(mat(X, modes), mat(Y, modes)) / sample.size, c(dimX, dimF)) + + # Compute standardized X and Y correlation tensor + SCovXY <- mlm(CovXY, Map(matpow, c(CovXX, CovYY), -1 / 2)) + + # mode-wise canonical correlation directions + hosvd <- HOSVD(SCovXY, nu = rep(pmin(dimX, dimF), 2L)) + dirsX <- hosvd$Us[modes] + dirsY <- hosvd$Us[modes + length(modes)] + + list(dirsX = dirsX, dirsY = dirsY) +} + + +# c(dirsX, dirsY) %<-% HOCCA(X, F) +# B.hocca <- Reduce(kronecker, rev(Map(tcrossprod, dirsX, dirsY))) +# dist.subspace(B.true, B.hocca, normalize = TRUE) +# dist.subspace(B.true, Reduce(kronecker, rev(dirsX)), normalize = TRUE) + +# cca <- cancor(mat(X, 4), mat(F, 4)) +# B.cca <- tcrossprod(cca$xcoef[, prod(dim(X)[-4])], cca$xcoef[, prod(dim(F)[-4])]) +# dist.subspace(B.true, cca$xcoef[, prod(dim(X)[-4])], normalize = TRUE) diff --git a/tensorPredictors/R/SIR.R b/tensorPredictors/R/SIR.R new file mode 100644 index 0000000..99cb02a --- /dev/null +++ b/tensorPredictors/R/SIR.R @@ -0,0 +1,36 @@ +#' Sliced Inverse Regression +#' +#' @export +SIR <- function(X, y, d, nr.slices = 10L, slice.method = c("cut", "ecdf")) { + + if (!(is.factor(y) || is.integer(y))) { + slice.method <- match.arg(slice.method) + if (slice.method == "ecdf") { + y <- cut(ecdf(y)(y), nr.slices) + } else { + y <- cut(y, nr.slices) + # ensure there are no empty slices + if (any(table(y) == 0)) { + y <- as.factor(as.integer(y)) + } + } + } + + # Center `X` + Z <- scale(X, scale = FALSE) + + # Split `Z` into slices determined by `y` + slices <- Map(function(i) Z[i, , drop = FALSE], split(seq_along(y), y)) + + # Sizes and Means for each slice + slice.sizes <- mapply(nrow, slices) + slice.means <- Map(colMeans, slices) + + # Inbetween slice covariances + sCov <- Reduce(`+`, Map(function(mean_s, n_s) { + n_s * tcrossprod(mean_s) + }, slice.means, slice.sizes)) / nrow(X) + + # Compute EDR directions + La.svd(sCov, d, 0L)$u +} diff --git a/tensorPredictors/R/det.R b/tensorPredictors/R/det.R new file mode 100644 index 0000000..215e4aa --- /dev/null +++ b/tensorPredictors/R/det.R @@ -0,0 +1,7 @@ +#' Determinant of a matrix +#' +#' @export +La.det <- function(A) { + storage.mode(A) <- "double" + .Call("C_det", A, PACKAGE = "tensorPredictors") +} diff --git a/tensorPredictors/R/gmlm_ising.R b/tensorPredictors/R/gmlm_ising.R index 108b47a..630fe21 100644 --- a/tensorPredictors/R/gmlm_ising.R +++ b/tensorPredictors/R/gmlm_ising.R @@ -4,7 +4,7 @@ #' #' @export gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), - proj.betas = NULL, proj.Omegas = NULL, + proj.betas = NULL, proj.Omegas = NULL, Omega.mask = NULL, max.iter = 1000L, eps = sqrt(.Machine$double.eps), step.size = 1e-3, @@ -112,7 +112,7 @@ gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), matX <- mat(X, sample.axis) degen <- crossprod(matX) == 0 degen.mask <- which(degen) - # If there are degenerate combination, compute an (arbitrary) bound the + # If there are degenerate combination, compute an (arbitrary) bound of the # log odds parameters of those combinations if (any(degen.mask)) { degen.ind <- arrayInd(degen.mask, dim(degen)) @@ -145,7 +145,7 @@ gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), } # Initialize mean squared gradients - grad2_betas <- Map(array, 0, Map(dim, betas)) + grad2_betas <- Map(array, 0, Map(dim, betas)) grad2_Omegas <- Map(array, 0, Map(dim, Omegas)) # Keep track of the last loss to accumulate loss difference sign changes @@ -166,6 +166,11 @@ gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), grad_betas <- Map(matrix, 0, dimX, dimF) Omega <- Reduce(kronecker, rev(Omegas)) + # Mask Omega, that is to enforce the "linear" constraint `T2` + if (!is.null(Omega.mask)) { + Omega[Omega.mask] <- 0 + } + # second order residuals accumulator # `sum_i (X_i o X_i - E[X o X | Y = y_i])` R2 <- array(0, dim = c(dimX, dimX)) @@ -186,7 +191,7 @@ gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), # accumulate loss matX_i <- mat(eval(`X[..., i]`), modes) loss <- loss - ( - sum(matX_i * (params_i %*% matX_i)) + n_i * log(attr(m2_i, "prob_0")) + sum(matX_i * (params_i %*% matX_i)) + n_i * attr(m2_i, "log_prob_0") ) R2_i <- tcrossprod(matX_i) - n_i * m2_i @@ -200,6 +205,12 @@ gmlm_ising <- function(X, F, y = NULL, sample.axis = length(dim(X)), R2 <- R2 + as.vector(R2_i) } + # Apply the `T2` constraint on the Residuals as well (refer to `T2`) + # That is, we compute G2 from g2 as in Theorem 2. + if (!is.null(Omega.mask)) { + R2[Omega.mask] <- 0 + } + grad_Omegas <- Map(function(j) { grad <- mlm(kronperm(R2), Map(as.vector, Omegas[-j]), modes[-j], transposed = TRUE) dim(grad) <- dim(Omegas[[j]]) diff --git a/tensorPredictors/R/ising_m2.R b/tensorPredictors/R/ising_m2.R index e85b5a4..b00160c 100644 --- a/tensorPredictors/R/ising_m2.R +++ b/tensorPredictors/R/ising_m2.R @@ -12,43 +12,7 @@ ising_m2 <- function( ) M2 <- vech.pinv(m2) - attr(M2, "prob_0") <- attr(m2, "prob_0") + attr(M2, "log_prob_0") <- attr(m2, "log_prob_0") M2 } - - -# library(tensorPredictors) - -# dimX <- c(3, 4) -# dimF <- rep(1L, length(dimX)) - -# betas <- Map(diag, 1, dimX, dimF) -# Omegas <- list( -# 1 - diag(dimX[1]), -# toeplitz(rev(seq(0, len = dimX[2])) / dimX[2]) -# ) -# Omega <- Reduce(kronecker, rev(Omegas)) - -# y <- array(1, dimF) -# params <- diag(as.vector(mlm(y, betas))) + Omega - -# # params <- array(0, dim(Omega)) - -# (prob_0 <- attr(ising_m2(params), "prob_0")) -# (probs <- replicate(20, attr(ising_m2(params, use_MC = TRUE, nr_threads = 8), "prob_0")))[1] -# m <- mean(probs) -# s <- sd(probs) - -# (prob_a <- (function(p, M2) { -# (1 + p * (p + 1) / 2 + 2 * sum(M2) - 2 * (p + 1) * sum(diag(M2))) / 2^p -# })(prod(dimX), ising_m2(params, use_MC = FALSE))) - -# par(mar = c(1, 2, 1, 2) + 0.1) -# plot(probs, ylim = pmax(0, range(probs, prob_0, prob_a)), pch = 16, cex = 1, -# xaxt = "n", xlab = "", col = "gray", log = "y", bty = "n") -# lines(cumsum(probs) / seq_along(probs), lty = 2, lwd = 2) -# abline(h = c(m - s, m, m + s), lty = c(3, 2, 3), col = "red", lwd = 2) -# abline(h = c(prob_0, prob_a), lwd = 2) -# axis(4, at = prob_0, labels = sprintf("%.1e", prob_0)) -# axis(4, at = prob_a, labels = sprintf("%.1e", prob_a)) diff --git a/tensorPredictors/R/kronperm.R b/tensorPredictors/R/kronperm.R new file mode 100644 index 0000000..4bf7fc6 --- /dev/null +++ b/tensorPredictors/R/kronperm.R @@ -0,0 +1,36 @@ +#' Kronecker Permutation of an array +#' +#' Computes a permutation and reshaping of `A` such that +#' kronperm(B %o% C) == kronecker(B, C) +#' +#' @param A multi-dimensional array +#' @param dims dimensions `A` should have overwriting the actuall dimensions +#' @param ncomp number of "components" counting the elements of an outer product +#' used to generate `A` if it is the result of an outer product. +#' +#' @examples +#' A <- array(rnorm(24), dim = c(2, 3, 4)) +#' B <- array(rnorm(15), dim = c(5, 3, 1)) +#' C <- array(rnorm(84), dim = c(7, 4, 3)) +#' +#' all.equal( +#' kronperm(outer(A, B)), +#' kronecker(A, B) +#' ) +#' all.equal( +#' kronperm(Reduce(outer, list(A, B, C)), ncomp = 3L), +#' Reduce(kronecker, list(A, B, C)) +#' ) +#' +#' @export +kronperm <- function(A, dims = dim(A), ncomp = 2L) { + # force `A` to have a multiple of `ncomp` dimensions + dim(A) <- c(dims, rep(1L, length(dims) %% ncomp)) + # compute axis permutation + perm <- as.vector(t(matrix(seq_along(dim(A)), ncol = ncomp)[, ncomp:1])) + # permute elements of A + K <- aperm(A, perm, resize = FALSE) + # collapse/set dimensions + dim(K) <- apply(matrix(dim(K), ncol = ncomp), 1, prod) + K +} diff --git a/tensorPredictors/R/mat_mani_projections.R b/tensorPredictors/R/mat_mani_projections.R new file mode 100644 index 0000000..b07138a --- /dev/null +++ b/tensorPredictors/R/mat_mani_projections.R @@ -0,0 +1,120 @@ +#' @rdname matProj +#' @export +projSym <- function(A) 0.5 * (A + t(A)) +#' @rdname matProj +#' @export +projDiag <- function(A) diag(diag(A)) +#' @rdname matProj +#' @export +.projBand <- function(dims, low, high) { + diag.index <- .row(dims) - .col(dims) + mask <- (diag.index <= low) & (-high <= diag.index) + function(A) A * mask +} +#' @rdname matProj +#' @export +.projSymBand <- function(dims, low, high) { + diag.index <- .row(dims) - .col(dims) + mask <- (diag.index <= low) & (-high <= diag.index) + function(A) projSym(A) * mask +} +#' @rdname matProj +#' @export +.projPSD <- function(sym = FALSE) { + if (sym) { + function(A) { + eig <- eigen(A, symmetric = TRUE) + eig$vectors %*% (pmax(0, eig$values) * t(eig$vectors)) + } + } else { + function(A) { + eig <- eigen(0.5 * (A + t(A)), symmetric = TRUE) + eig$vectors %*% (pmax(0, eig$values) * t(eig$vectors)) + } + } +} +#' @rdname matProj +#' @export +.projRank <- function(rank) { + force(rank) + function(A) { + rank <- min(dim(A), rank) + svdA <- La.svd(A, rank, rank) + svdA$u %*% (svdA$d[seq_len(rank)] * svdA$vt) + } +} +#' @rdname matProj +#' @export +.projSymRank <- function(rank) { + force(rank) + function(A) { + rank <- min(dim(A), rank) + svdA <- La.svd(0.5 * (A + t(A)), rank, rank) + svdA$u %*% (svdA$d[seq_len(rank)] * svdA$vt) + } +} +#' @rdname matProj +#' @export +projStiefel <- function(A) { + # Using a polar decomposition of `A = Q P` via SVD `A = U D V^T`. Compaired + # to a QR decomposition the polar decomposition is unique, making it "stabel". + svdA <- La.svd(A) + svdA$u %*% svdA$vt # = Q +} +# .projKron <- function(dims) { +# ... # TODO: Implement this! +# } + +#' @rdname matProj +#' @export +.projMaskedMean <- function(mask) { + force(mask) + function(A) { + `[<-`(matrix(0, nrow(A), ncol(A)), mask, mean(A[mask])) + } +} + +#' Projections onto matrix manifolds +#' +#' @examples +#' p <- 5 +#' q <- 4 +#' A <- matrix(rnorm(p * q), p, q) +#' +#' # General Matrices +#' matProj("TriDiag", dim(A))(A) +#' matProj("Band", dim(A), low = 1, high = 2)(A) +#' matProj("Rank", rank = 2)(A) +#' matProj("Stiefel")(A) +#' +#' # Symmetric projections need square matrices +#' S <- matrix(rnorm(p^2), p) +#' +#' matProj("Sym")(S) +#' matProj("SymTriDiag", dim(S))(S) +#' matProj("SymBand", dim(S), low = 1, high = 2)(S) +#' matProj("PSD")(S) +#' matProj("SymRank", rank = 1)(S) +#' +#' @rdname matProj +#' +#' @export +matProj <- function(manifold, dims = NULL, low = NULL, high = NULL, sym = FALSE, rank = NULL) { + + switch(tolower(manifold), + identity = identity, + sym = projSym, + tridiag = .projBand(dims, 1L, 1L), + symtridiag = .projSymBand(dims, 1L, 1L), + band = .projBand(dims, low, high), + symband = .projSymBand(dims, low, high), + psd = .projPSD(sym), + rank = .projRank(rank), + symrank = .projSymRank(rank), + stiefel = projStiefel + ) + +} + +# #' Basis of .... +# mat.proj.basis <- function(manifold, dims = NULL, low = NULL, high = NULL, sym = FALSE, rank = NULL) ... diff --git a/tensorPredictors/R/pinv.R b/tensorPredictors/R/pinv.R new file mode 100644 index 0000000..54f2649 --- /dev/null +++ b/tensorPredictors/R/pinv.R @@ -0,0 +1,17 @@ +#' Moore-Penrose Pseudo inverse +#' +#' @param A any matrix +#' +#' @returns another matrix +#' +#' @export +pinv <- function(A) { + A <- as.matrix(A) + if (nrow(A) < ncol(A)) { + crossprod(A, matpow(tcrossprod(A), -1)) + } else if (nrow(A) > ncol(A)) { + tcrossprod(matpow(crossprod(A), -1), A) + } else { + matpow(A, -1) + } +} diff --git a/tensorPredictors/R/slice_select.R b/tensorPredictors/R/slice_select.R new file mode 100644 index 0000000..645809a --- /dev/null +++ b/tensorPredictors/R/slice_select.R @@ -0,0 +1,67 @@ +#' Slice index selection +#' +#' @examples +#' # Exquivalent to +#' array(A[slice.index(A, mode) == index], dim = dim(A)[-mode]) +#' +#' @export +slice.select <- function(A, mode, index) { + arg <- rep("", length(dim(A))) + arg[mode] <- "i" + expr <- str2lang(paste0("A[", paste0(arg, collapse = ","), "]", collapse = "")) + slice <- eval(expr, list(i = index)) + dim(slice) <- dim(A)[-mode] + slice +} + +#' +#' @export +slice.expr <- function(A, mode, index = "i", drop = TRUE, nr.axis = length(dim(A))) { + str <- as.character(substitute(A)) + arg <- rep("", nr.axis) + arg[mode] <- as.character(substitute(index)) + str2lang(paste0(str, "[", paste0(arg, collapse = ","), + if (drop) "]" else ",drop=FALSE]", collapse = "")) +} + + +#' @export +slice.assign.expr <- function(obj, nr.axis) { + assign.call <- as.call(c( + list(`[<-`, substitute(obj)), + rep(list(alist(a = )$a), nr.axis - 1L), # replicate empty symbol + substitute(index), substitute(x) + )) + function(i, val) { + eval(assign.call, envir = list(index = i, x = val)) + } +} + + +# n <- 1000 +# p <- c(2, 4, 3) +# A <- array(seq_len(prod(n, p)), dim = c(p, n)) + +# mode <- 4 +# index <- 7 + +# stopifnot(all.equal( +# A[, , , index], +# array(A[slice.index(A, mode) == index], dim = dim(A)[-mode]) +# )) +# stopifnot(all.equal( +# A[, , , index], +# slice.select(A, mode, index) +# )) + +# arg <- rep("", length(dim(A))) +# arg[mode] <- "i" +# `A[..., i]` <- str2lang(paste0("A[", paste0(arg, collapse = ","), "]", collapse = "")) + +# microbenchmark::microbenchmark( +# A[, , , index], +# eval(`A[..., i]`, list(i = index)), +# slice.select(A, mode, index), +# array(A[slice.index(A, mode) == index], dim = dim(A)[-mode]) +# ) + diff --git a/tensorPredictors/R/solve.R b/tensorPredictors/R/solve.R new file mode 100644 index 0000000..f718f20 --- /dev/null +++ b/tensorPredictors/R/solve.R @@ -0,0 +1,11 @@ +#' Linear Equation Solver +#' +#' Using Lapack DGESV, similar to the base solve routine of R. +#' +#' @note for testing purposes +#' +#' @export +La.solve <- function(A, B = diag(nrow(A))) { + storage.mode(A) <- storage.mode(B) <- "double" + .Call("C_solve", A, as.matrix(B), PACKAGE = "tensorPredictors") +} diff --git a/tensorPredictors/R/sym.R b/tensorPredictors/R/sym.R new file mode 100644 index 0000000..752e8d3 --- /dev/null +++ b/tensorPredictors/R/sym.R @@ -0,0 +1,81 @@ +#' Generale a matrix of all permutations of `n` elements +permutations <- function(n) { + if (n <= 0) { + matrix(nrow = 0, ncol = 0) + } else if (n == 1) { + matrix(1) + } else { + sub.perm <- permutations(n - 1) + p <- nrow(sub.perm) + A <- matrix(NA, n * p, n) + for (i in 1:n) { + A[(i - 1) * p + 1:p, ] <- cbind(i, sub.perm + (sub.perm >= i)) + } + A + } +} + + +#' General symmetrization opperation for tensors (arrays) of equal dimensions +#' +#' @param A array of dimensions c(p, ..., p) +#' +#' @returns array of same dimensions as `A` +#' +#' @export +tsym <- function(A) { + stopifnot(all(dim(A) == nrow(A))) + + if (is.matrix(A)) { + return(0.5 * (A + t(A))) + } + + axis.perm <- permutations(length(dim(A))) + + S <- array(0, dim(A)) + for (i in seq_len(nrow(axis.perm))) { + S <- S + aperm(A, axis.perm[i, ]) + } + + S / nrow(axis.perm) +} + + + +#' Genralized (pseudo) symmetrication for generel multi-dimensional arrays +sym <- function(A, FUN = `+`, scale = factorial(length(dim(A)))) { + FUN <- match.fun(FUN) + + if (is.matrix(A) && (nrow(A) == ncol(A))) { + A <- FUN(A, t(A)) + return(if (is.numeric(scale)) A / scale else A) + } + + A.copy <- A + perm <- seq_along(dim(A)) + while (length(pivot <- which(diff(perm) > 0))) { + pivot <- max(pivot) + successor <- max(which(perm[seq_along(perm) > pivot] > perm[pivot])) + pivot + perm[c(pivot, successor)] <- perm[c(successor, pivot)] + suffix <- seq(pivot + 1, length(perm)) + perm <- c(perm[-suffix], perm[rev(suffix)]) + + modes <- which(perm != seq_along(perm)) + sub.dimA <- dim(A) + sub.dimA[modes] <- min(dim(A)[modes]) + sub.indices <- Map(seq_len, sub.dimA) + sub.selection <- do.call(`[`, c(list(A.copy), sub.indices, drop = FALSE)) + + + sub.assign <- do.call(call, c(list("[<-", quote(quote(A))), sub.indices, + quote(quote(FUN( + do.call(`[`, c(list(A), sub.indices)), + aperm(do.call(`[`, c(list(A.copy), sub.indices, drop = FALSE)), perm) + ))) + )) + + A <- eval(sub.assign) + } + + if (is.numeric(scale)) A / scale else A +} diff --git a/tensorPredictors/src/det.c b/tensorPredictors/src/det.c new file mode 100644 index 0000000..b355e0a --- /dev/null +++ b/tensorPredictors/src/det.c @@ -0,0 +1,73 @@ +#include "det.h" + + +// For reference see: `det_ge_real` in "R-4.2.1/src/modules/lapack/Lapack.c" +// of the R source code. +double det( + /* dim */ const int dimA, + /* matrix */ const double* A, const int ldA, + double* work_mem, int* info +) { + // if working memory size query, return immediately + if (work_mem == NULL) { + *info = dimA * (dimA + 1); + return 0.0; + } + + // determinant of "zero size" matrix is 1 (by definition) + if (dimA == 0) { + return 1.0; + } + + // copy `A` to (continuous) `work_mem` cause `dgetrf` works "in place" + for (int i = 0; i < dimA; ++i) { + memcpy(work_mem + i * dimA, A + i * ldA, dimA * sizeof(double)); + } + + // L U factorization of `A` + int error = 0; + int* ipvt = (int*)(work_mem + dimA * dimA); + F77_CALL(dgetrf)(&dimA, &dimA, work_mem, &dimA, ipvt, &error); + + // check if an error occured which is the case iff `dgetrf` gives a negative error + *info |= (error < 0) * error; + // in both cases we return zero (ether the determinant is zero or error occured) + if (error) { + return 0.0; + } + + // res <- det(A) = sign(P) * prod(diag(L)) where P is the pivoting permutation + double res = 1.0; + for (int i = 0; i < dimA; ++i) { + res *= (ipvt[i] != (i + 1) ? -1.0 : 1.0) * work_mem[i * (dimA + 1)]; + } + + return res; +} + +/** + * R bindong to `det` + */ +extern SEXP R_det(SEXP A) { + // check if A is a real valued square matrix + if (!Rf_isReal(A) || !Rf_isMatrix(A) || Rf_nrows(A) != Rf_ncols(A)) { + Rf_error("`A` must be a real valued squae matrix"); + } + + // allocate working memory + int work_size; + (void)det(Rf_nrows(A), NULL, Rf_nrows(A), NULL, &work_size); + double* work_mem = (double*)R_alloc(work_size, sizeof(double)); + + // compute determinant (followed only by if-statement, no protection required) + int error = 0; + SEXP res = Rf_ScalarReal( + det(Rf_nrows(A), REAL(A), Rf_nrows(A), work_mem, &error) + ); + + if (error) { + Rf_error("Encountered error code %d in `det`", error); + } + + return res; +} diff --git a/tensorPredictors/src/det.h b/tensorPredictors/src/det.h new file mode 100644 index 0000000..25ba156 --- /dev/null +++ b/tensorPredictors/src/det.h @@ -0,0 +1,15 @@ +#ifndef INCLUDE_GUARD_DET_H +#define INCLUDE_GUARD_DET_H + +#include "R_api.h" + +/** + * Determinant of a matrix (or log of determinant) + */ +double det( + /* dim */ const int dimA, + /* matrix */ const double* A, const int ldA, + double* work_mem, int* info +); + +#endif /* INCLUDE_GUARD_DET_H */ diff --git a/tensorPredictors/src/ising_m2.c b/tensorPredictors/src/ising_m2.c index dbf0f8a..6d7f02a 100644 --- a/tensorPredictors/src/ising_m2.c +++ b/tensorPredictors/src/ising_m2.c @@ -27,6 +27,8 @@ * [ * * * * * * ] [ * * * 17 * * ] [ * ] * */ +#include // DBL_MAX + #include "R_api.h" #include "bit_utils.h" #include "int_utils.h" @@ -121,7 +123,7 @@ double ising_m2_exact(const size_t dim, const double* params, double* M2) { const double prob_X = exp(dot_X); sum_0 += prob_X; - // Accumulate set bits probability for the first end second moment `E[X X']` + // Accumulate set bits probability for the first and second moment `E[X X']` for (uint32_t Y = X; Y; Y &= Y - 1) { const int i = bitScanLS32(Y); const int I = (i * (2 * dim - 1 - i)) / 2; @@ -137,7 +139,7 @@ double ising_m2_exact(const size_t dim, const double* params, double* M2) { M2[i] *= prob_0; } - return prob_0; + return log(prob_0); } @@ -162,7 +164,7 @@ double ising_m2_MC( int* X = (int*)R_alloc(dim, sizeof(int)); // Accumulator for Monte-Carlo estimate for zero probability `P(X = 0)` - double accum = 0.0; + double accum = 0.0, max_mdot_X = -DBL_MAX; // Create/Update R's internal PRGN state GetRNGstate(); @@ -182,7 +184,11 @@ double ising_m2_MC( dot_X += (X[i] & X[j]) ? params[I + j] : 0.0; } } - accum += exp(-dot_X); + if (-dot_X > max_mdot_X) { + accum *= exp(max_mdot_X + dot_X); + max_mdot_X = -dot_X; + } + accum += exp(-dot_X - max_mdot_X); } // Write R's internal PRNG state back (if needed) @@ -190,11 +196,11 @@ double ising_m2_MC( // Compute means from counts for (size_t i = 0; i < len; ++i) { - M2[i] = (double)counts[i] / (double)(nr_samples); + M2[i] = (double)counts[i] / (double)nr_samples; } // Prob. of zero event (Ising p.m.f. scaling constant) - return accum / (exp2((double)dim) * (double)nr_samples); + return log(accum) + max_mdot_X - log(2) * (double)dim - log((double)nr_samples); } @@ -210,7 +216,7 @@ typedef struct thrd_data { const double* params; // Ising model parameters int* X; // Working memory to store current binary sample uint32_t* counts; // (output) count of single and two way interactions - double accum; // (output) Monte-Carlo accumulator for `P(X = 0)` + double log_sum_exp; // (output) Monte-Carlo inbetween value for `log(P(X = 0))` } thrd_data_t; // Worker thread function @@ -222,8 +228,8 @@ int thrd_worker(thrd_data_t* data) { // Initialize counts to zero (void)memset(data->counts, 0, (dim * (dim + 1) / 2) * sizeof(uint32_t)); - // Init Monte-Carlo estimate for zero probability `P(X = 0)` - data->accum = 0.0; + // Accumulator for Monte-Carlo estimate for zero probability `P(X = 0)` + double accum = 0.0, max_mdot_X = -DBL_MAX; // Spawn Monte-Carlo Chains (one foe every sample) for (size_t sample = 0; sample < data->nr_samples; ++sample) { @@ -240,8 +246,14 @@ int thrd_worker(thrd_data_t* data) { dot_X += (X[i] & X[j]) ? data->params[I + j] : 0.0; } } - data->accum += exp(-dot_X); - } + if (-dot_X > max_mdot_X) { + accum *= exp(max_mdot_X + dot_X); + max_mdot_X = -dot_X; + } + accum += exp(-dot_X - max_mdot_X); + } + + data->log_sum_exp = log(accum) + max_mdot_X; return 0; } @@ -286,13 +298,19 @@ double ising_m2_MC_thrd( } // Accumulate worker results into first (tid = 0) worker counts result - // as well as accumulate the Monte-Carlo accumulators into one - double accum = threads_data[0].accum; + // while determining the log sum exp maximum + double max = threads_data[0].log_sum_exp; for (size_t tid = 1; tid < nr_threads; ++tid) { for (size_t i = 0; i < len; ++i) { counts[i] += counts[tid * len + i]; } - accum += threads_data[tid].accum; + max = (max < threads_data[tid].log_sum_exp) ? threads_data[tid].log_sum_exp : max; + } + + // Accum all `log(P(X = 0))` via "LogSumExp" trick into final result + double accum = 0; + for (size_t tid = 0; tid < nr_threads; ++tid) { + accum += exp(threads_data[tid].log_sum_exp - max); } // convert discreat counts into means @@ -300,8 +318,8 @@ double ising_m2_MC_thrd( M2[i] = (double)counts[i] / (double)nr_samples; } - // Prob. of zero event (Ising p.m.f. scaling constant) - return accum / (exp2((double)dim) * (double)nr_samples); + // Log of Prob. of zero event (Ising p.m.f. scaling constant) + return log(accum) + max - log(2) * (double)dim - log((double)nr_samples); } #endif /* !__STDC_NO_THREADS__ */ @@ -363,9 +381,9 @@ extern SEXP R_ising_m2( SEXP _M2 = PROTECT(Rf_allocVector(REALSXP, dim * (dim + 1) / 2)); ++protect_count; - // asside computed zero event probability (inverse partition function), the - // scaling factor for the Ising model p.m.f. - double prob_0 = -1.0; + // asside computed log of zero event probability (log of recibrocal partition + // function), the log scaling factor for the Ising model p.m.f. + double log_prob_0 = -1.0; if (use_MC) { // Convert and validate arguments for the Monte-Carlo methods @@ -384,15 +402,15 @@ extern SEXP R_ising_m2( if (nr_threads == 1) { // Single threaded Monte-Carlo method - prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); + log_prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); } else { // Multi-Threaded Monte-Carlo method if provided, otherwise use // the single threaded version with a warning #ifdef __STDC_NO_THREADS__ Rf_warning("Multi-Threading NOT supported, using fallback."); - prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); + log_prob_0 = ising_m2_MC(nr_samples, warmup, dim, params, REAL(_M2)); #else - prob_0 = ising_m2_MC_thrd( + log_prob_0 = ising_m2_MC_thrd( nr_samples, warmup, nr_threads, dim, params, REAL(_M2) ); @@ -405,13 +423,13 @@ extern SEXP R_ising_m2( } // and call the exact method - prob_0 = ising_m2_exact(dim, params, REAL(_M2)); + log_prob_0 = ising_m2_exact(dim, params, REAL(_M2)); } // Set log-lokelihood as an attribute to the computed second moment - SEXP _prob_0 = PROTECT(Rf_ScalarReal(prob_0)); + SEXP _log_prob_0 = PROTECT(Rf_ScalarReal(log_prob_0)); ++protect_count; - Rf_setAttrib(_M2, Rf_install("prob_0"), _prob_0); + Rf_setAttrib(_M2, Rf_install("log_prob_0"), _log_prob_0); // release SEPXs to the garbage collector UNPROTECT(protect_count); diff --git a/tensorPredictors/src/solve.c b/tensorPredictors/src/solve.c new file mode 100644 index 0000000..9a3405b --- /dev/null +++ b/tensorPredictors/src/solve.c @@ -0,0 +1,100 @@ +#include "solve.h" + +void solve( + /* dims */ const int dimA, const int nrhs, + /* matrix */ const double* A, const int ldA, + /* matrix */ const double* B, const int ldB, + /* matrix */ double* X, const int ldX, + double* work_mem, int* info +) { + // Compute required working memory size if requested + if (work_mem == NULL) { + *info = dimA * (dimA + 1); + return; + } + + // Copy `A` to (continuous) working memory + for (int i = 0; i < dimA; ++i) { + memcpy(work_mem + i * dimA, A + i * ldA, dimA * sizeof(double)); + } + + // Copy `B` to `X` or set `X` to identity + if (B == NULL) { + double* X_col = X; + for (int j = 0; j < dimA; ++j, X_col += ldX) { + for (int i = 0; i < dimA; ++i) { + *(X_col + i) = (i == j) ? 1.0 : 0.0; + } + } + } else { + for (int i = 0; i < nrhs; ++i) { + memcpy(X + i * ldX, B + i * ldB, dimA * sizeof(double)); + } + } + + // Lapack routine DGESV to solve linear system A X = B which writes + // result into `A`, `B` which are copied into working memory and the result + // memory `X` + int error = 0; + F77_CALL(dgesv)( + /* dims */ &dimA, &nrhs, + /* matrix A */ work_mem, &dimA, /* [in,out] A -> P L U */ + /* ipiv */ (int*)(work_mem + dimA * dimA), /* [out] */ + /* matrix B */ X, &ldX, /* [in,out] B -> X */ + &error /* [out] */ + ); + + // update error flag + *info |= error; +} + +/** + * R binding to `solve` which solves A X = B for X + */ +extern SEXP R_solve(SEXP A, SEXP B) { + // Check types + if (!(Rf_isReal(A) && Rf_isMatrix(A)) + || !(Rf_isReal(B) && Rf_isMatrix(B))) { + Rf_error("All arguments must be real valued matrices"); + } + + // check dimensions + if (Rf_nrows(A) != Rf_ncols(A) + || Rf_ncols(A) != Rf_nrows(B)) { + Rf_error("Dimension missmatch"); + } + + // Allocate result matrix `X` + SEXP X = PROTECT(Rf_allocMatrix(REALSXP, Rf_nrows(B), Rf_ncols(B))); + + // Allocate required working memory + int work_size = 0; + solve( + Rf_nrows(A), Rf_ncols(B), + NULL, Rf_nrows(A), + NULL, Rf_nrows(B), + NULL, Rf_nrows(X), + NULL, &work_size + ); + double* work_mem = (double*)R_alloc(work_size, sizeof(double)); + + // Solve the system A X = B an write results into `X` + int error = 0; + solve( + Rf_nrows(A), Rf_ncols(B), + REAL(A), Rf_nrows(A), + REAL(B), Rf_nrows(B), + REAL(X), Rf_nrows(X), + work_mem, &error + ); + + // release `X` to the garbage collector + UNPROTECT(1); + + // check error after unprotect + if (error) { + Rf_error("Solve ended with error code %d", error); + } + + return X; +} diff --git a/tensorPredictors/src/solve.h b/tensorPredictors/src/solve.h new file mode 100644 index 0000000..e8e623e --- /dev/null +++ b/tensorPredictors/src/solve.h @@ -0,0 +1,19 @@ +#ifndef INCLUDE_GUARD_SOLVE_H +#define INCLUDE_GUARD_SOLVE_H + +#include "R_api.h" + +/** + * Solves a linear equation system of the form + * A X = B + * for `X` where `A` is a `dimA` x `dimA` matrix and `B` is `dimA` x `nrhs`. + */ +void solve( + /* dims */ const int dimA, const int nrhs, + /* matrix */ const double* A, const int ldA, + /* matrix */ const double* B, const int ldB, + /* matrix */ double* X, const int ldX, + double* work_mem, int* info +); + +#endif /* INCLUDE_GUARD_SOLVE_H */