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

@@ -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")

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)
+}

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)

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
+}

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")
+}

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))
@@ -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]])

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))

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
+}

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) ...

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)
+    }
+}

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])
+# )
+

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")
+}

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
+}

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;
+}

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 */

tensorPredictors/src/ising_m2.c

@@ -27,6 +27,8 @@
  * [ *    *    *    *    *    * ]    [ *    *    *   17    *    * ]   [ * ]
  *
  */
+#include <float.h>  // 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)`
+    // Accumulator for Monte-Carlo estimate for zero probability `P(X = 0)`
-    data->accum = 0.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
+    // while determining the log sum exp maximum
-    double accum = threads_data[0].accum;
+    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)
+    // Log of Prob. of zero event (Ising p.m.f. scaling constant)
-    return accum / (exp2((double)dim) * (double)nr_samples);
+    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
+    // asside computed log of zero event probability (log of recibrocal partition
-    // scaling factor for the Ising model p.m.f.
+    // function), the log scaling factor for the Ising model p.m.f.
-    double prob_0 = -1.0;
+    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);

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;
+}

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 */