## -----------------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

## -----------------------------------------------------------------------------
library(neuromapr)

## -----------------------------------------------------------------------------
set.seed(1)
neuromapr:::random_rotation_euler()

## -----------------------------------------------------------------------------
set.seed(1)
neuromapr:::random_rotation_rodrigues()

## -----------------------------------------------------------------------------
set.seed(42)
n_lh <- 50
n_rh <- 50
coords <- list(
  lh = matrix(rnorm(n_lh * 3), ncol = 3),
  rh = matrix(rnorm(n_rh * 3), ncol = 3)
)
coords$lh <- t(apply(coords$lh, 1, function(x) x / sqrt(sum(x^2))))
coords$rh <- t(apply(coords$rh, 1, function(x) x / sqrt(sum(x^2))))

data_x <- rnorm(n_lh + n_rh)
data_y <- 0.3 * data_x + rnorm(n_lh + n_rh, sd = 0.9)

## -----------------------------------------------------------------------------
nulls_euler <- null_alexander_bloch(
  data_x,
  coords,
  n_perm = 500L,
  seed = 1,
  rotation = "euler"
)

nulls_rodrigues <- null_alexander_bloch(
  data_x,
  coords,
  n_perm = 500L,
  seed = 1,
  rotation = "rodrigues"
)

## -----------------------------------------------------------------------------
null_cors_euler <- apply(nulls_euler$nulls, 2, cor, data_y)
null_cors_rodrigues <- apply(nulls_rodrigues$nulls, 2, cor, data_y)

df <- data.frame(
  r = c(null_cors_euler, null_cors_rodrigues),
  method = rep(c("Euler (ZYZ)", "Rodrigues (axis-angle)"), each = 500)
)

ggplot2::ggplot(df, ggplot2::aes(x = r, fill = method)) +
  ggplot2::geom_density(alpha = 0.5) +
  ggplot2::scale_fill_manual(values = c("steelblue", "darkorange")) +
  ggplot2::labs(x = "Null correlation (r)", y = "Density", fill = "Method") +
  ggplot2::theme_minimal()

## -----------------------------------------------------------------------------
obs_r <- cor(data_x, data_y)
p_euler <- mean(abs(null_cors_euler) >= abs(obs_r))
p_rodrigues <- mean(abs(null_cors_rodrigues) >= abs(obs_r))

data.frame(
  method = c("euler", "rodrigues"),
  p_value = c(p_euler, p_rodrigues)
)

## -----------------------------------------------------------------------------
# null_alexander_bloch(
#   data,
#   coords,
#   n_perm = 1000L,
#   rotation = "rodrigues"
# )