## ----setup, include = FALSE---------------------------------------------------
options(rmarkdown.html_vignette.check_title = FALSE)
knitr::opts_chunk$set(
  collapse = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  comment = "#>"
)

# Optional packages required by the regression and ridge chunks
has_sc3 <- requireNamespace("scatterplot3d", quietly = TRUE)
has_sdetorus <- requireNamespace("sdetorus", quietly = TRUE)
has_smacof <- requireNamespace("smacof", quietly = TRUE)

## ----load---------------------------------------------------------------------
library(polykde)

## ----data-layout--------------------------------------------------------------
d <- c(1, 2)
r <- length(d)
X <- r_unif_polysph(n = 4, d = d)
X
dim(X)

# The 0-based indexes delimiting the blocks: columns 1:2 and 3:5
comp_ind_dj(d = d)

# Each block has unit norm
rowSums(X[, 1:2]^2)
rowSums(X[, 3:5]^2)

## ----data-angles--------------------------------------------------------------
theta <- rbind(c(pi / 2, pi / 4), c(pi / 3, -pi / 2))
x <- angles_to_sph(theta = theta)
x
sph_to_angles(x = x)

## ----sim-vmf------------------------------------------------------------------
set.seed(42)
d <- c(2, 2)
mu <- c(0, 0, 1, 1, 0, 0)
kappa <- c(5, 10)
X <- r_vmf_polysph(n = 500, d = d, mu = mu, kappa = kappa)

# Densities at the modal point
d_vmf_polysph(x = rbind(mu), d = d, mu = mu, kappa = kappa)
d_unif_polysph(x = rbind(mu), d = d)

## ----sim-mixture--------------------------------------------------------------
mu_mix <- rbind(c(0, 0, 1, 1, 0, 0), c(1, 0, 0, 0, 0, 1))
kappa_mix <- rbind(c(10, 10), c(20, 5))
prop <- c(0.7, 0.3)
X_mix <- r_mvmf_polysph(n = 500, d = d, mu = mu_mix, kappa = kappa_mix,
                        prop = prop)
head(d_mvmf_polysph(x = X_mix, d = d, mu = mu_mix, kappa = kappa_mix,
                    prop = prop))

## ----sim-path-----------------------------------------------------------------
Y <- r_path_s2r(n = 200, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
dim(Y)

## ----kern-plot, fig.asp = 0.7, fig.width = 6----------------------------------
t_grid <- seq(0, 2, l = 200)
plot(t_grid, L(t = t_grid, kernel = "1"), type = "l", ylim = c(0, 1), lwd = 2,
     xlab = "t", ylab = "L(t)")
lines(t_grid, L(t = t_grid, kernel = "2"), col = 2, lwd = 2)
lines(t_grid, L(t = t_grid, kernel = "3", k = 10), col = 3, lwd = 2)
legend("topright", legend = c("von Mises-Fisher", "Epanechnikov", "softplus"),
       lwd = 2, col = 1:3)

# Normalizing constant of the product vMF kernel on (S^2)^2
c_kern(h = c(0.2, 0.2), d = c(2, 2), kernel = 1, kernel_type = 1)

## ----kern-eff-----------------------------------------------------------------
eff_kern(d = 2, r = 2, kernel = "1", kernel_type = "sph")
eff_kern(d = 2, r = 2, kernel = "3", k = 10, kernel_type = "sph")

## ----kde-sphere, fig.asp = 0.55, fig.width = 7--------------------------------
set.seed(42)
d <- 2
X <- r_mvmf_polysph(n = 500, d = d, mu = rbind(c(0, 0, 1),
                                               c(1, 1, 0) / sqrt(2)),
                    kappa = rbind(10, 25), prop = c(0.6, 0.4))
h <- bw_rot_polysph(X = X, d = d)$bw

# Evaluation on a regular grid of spherical angles, with theta_1 the polar
# angle and theta_2 the azimuth (see sph_to_angles())
th_1 <- seq(0, pi, l = 100)
th_2 <- seq(0, 2 * pi, l = 200)
grid <- angles_to_sph(theta = as.matrix(expand.grid(th_1, th_2)))
kde <- matrix(kde_polysph(x = grid, X = X, d = d, h = h), nrow = length(th_1))
image(th_1, th_2, kde, col = hcl.colors(50, "Viridis"), useRaster = TRUE,
      xlab = expression(theta[1]), ylab = expression(theta[2]))
contour(th_1, th_2, kde, add = TRUE, col = "white", drawlabels = FALSE)
points(sph_to_angles(x = X), pch = 16, cex = 0.4)

## ----kde-cv-------------------------------------------------------------------
head(r_kde_polysph(n = 3, X = X, d = d, h = h))
log_dens <- log_cv_kde_polysph(X = X, d = d, h = h)
summary(drop(log_dens))

## ----bw-compare---------------------------------------------------------------
set.seed(423432)
d <- c(1, 2)
mu <- c(0, 1, 0, 0, 1)
kappa <- c(5, 10)
X <- r_vmf_polysph(n = 200, d = d, mu = mu, kappa = kappa)

bw_mrot_polysph(X = X, d = d)
bw_rot_polysph(X = X, d = d)$bw
bw_cv_polysph(X = X, d = d, type = "LCV")$bw

## ----bw-lscv------------------------------------------------------------------
bw_cv_polysph(X = X, d = d, type = "LSCV", exact_vmf = TRUE)$bw
bw_cv_polysph(X = X, d = d, type = "LSCV", exact_vmf = TRUE,
              common_h = TRUE)$bw

## ----kre, fig.asp = 1, fig.width = 5, eval = has_sc3 && has_sdetorus----------
set.seed(423432)
n <- 100
x <- seq(0, 1, l = n)
Y <- r_path_s2r(n = n, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
h_kre <- bw_cv_kre_polysph(X = x, Y = Y, d = 2, p = 0, plot_cv = FALSE)$h_1se
h_kre

x_grid <- seq(-0.25, 1.25, l = 200)
sc3 <- scatterplot3d::scatterplot3d(Y, pch = 16, xlim = c(-1, 1),
                                    ylim = c(-1, 1), zlim = c(-1, 1),
                                    xlab = "", ylab = "", zlab = "",
                                    color = gray(0.5))
sc3$points3d(kre_polysph(x = x_grid, X = x, Y = Y, d = 2, h = h_kre, p = 0),
             type = "l", col = 2, lwd = 2)
sc3$points3d(kre_polysph(x = x_grid, X = x, Y = Y, d = 2, h = h_kre, p = 1),
             type = "l", col = 3, lwd = 2, lty = 2)
legend("topleft", legend = c("Local constant", "Local linear"), lwd = 2,
       lty = 1:2, col = 2:3)

## ----ridge, fig.asp = 1, fig.width = 5, eval = has_sc3 && has_sdetorus--------
set.seed(987204452)
n <- 100
X <- r_path_s2r(n = n, r = 1, spiral = FALSE, Theta = cbind(c(1, 0, 0)),
                sigma = 0.35)[, , 1]
rid <- euler_ridge(x = X, X = X, d = 2, h = 0.5, h_euler = 0.25, N = 30,
                   keep_paths = TRUE, show_prog = FALSE)

# The ridge points, the iterations spent, and the convergence flags
head(rid$ridge_y)
table(rid$iter)
table(rid$conv)

# Starting points, ridge points, and the paths joining them
sc3 <- scatterplot3d::scatterplot3d(X, pch = 16, xlim = c(-1, 1),
                                    ylim = c(-1, 1), zlim = c(-1, 1),
                                    xlab = "x", ylab = "y", zlab = "z",
                                    color = gray(0.5))
for (k in seq_len(nrow(rid$paths))) {

  sc3$points3d(t(rid$paths[k, , ]), type = "l", col = gray(0.75))

}
sc3$points3d(rid$ridge_y, pch = 16, cex = 0.75, col = 2)

## ----ridge-index, eval = has_sc3 && has_sdetorus && has_smacof----------------
rid_clean <- clean_euler_ridge(e = rid, X = X)
ind <- index_ridge(endpoints = rid_clean$ridge_y, X = X, d = 2,
                   probs_scores = seq(0, 1, l = 50))
summary(ind$scores_X)

## ----hom-test-----------------------------------------------------------------
set.seed(423432)
n <- c(100, 100)
X1 <- r_vmf_polysph(n = n[1], d = 2, mu = c(0, 0, 1), kappa = 5)
X2 <- r_vmf_polysph(n = n[2], d = 2, mu = c(0, 0, 1), kappa = 5)

# Homogeneity holds
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
                 type = "jsd", h = 0.5, B = 500, show_prog = FALSE)

# It does not: the second sample is more concentrated
X2 <- r_vmf_polysph(n = n[2], d = 2, mu = c(0, 0, 1), kappa = 20)
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
                 type = "jsd", h = 0.5, B = 500, show_prog = FALSE)

# The mean-based test is blind to this alternative
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
                 type = "mean", B = 500, show_prog = FALSE)

## ----hippo-data---------------------------------------------------------------
data("hippocampus")
dirs <- hippocampus$dirs
r <- ncol(dirs)
d <- rep(2, r)
X <- do.call(cbind, args = lapply(1:r, function(i) dirs[, i, ]))
dim(X)

## ----hippo-depth--------------------------------------------------------------
h_mrot <- bw_mrot_polysph(X = X, d = d, kernel = 3, k = 100)
h_rot <- bw_rot_polysph(X = X, d = d, bw0 = 2^(-1:2) %o% h_mrot,
                        kernel = 3, kernel_type = 2, k = 100)$bw
log_dens <- log_cv_kde_polysph(X = X, d = d, h = h_rot, wrt_unif = TRUE,
                               kernel = 3, kernel_type = 2, k = 100)
deep <- which.max(log_dens)
outly <- which.min(log_dens)
c("deepest" = deep, "most outlying" = outly)

## ----hippo-srep, fig.asp = 0.5, fig.width = 7, eval = requireNamespace("scatterplot3d", quietly = TRUE)----
par(mfrow = c(1, 2))
titles <- c("Deepest", "Most outlying")
for (j in 1:2) {

  i <- c(deep, outly)[j]
  view_srep(base = hippocampus$base[i, , ], dirs = dirs[i, , ],
            bdry = hippocampus$bdry[i, , ], static = TRUE, main = titles[j])

}

## ----hippo-test---------------------------------------------------------------
set.seed(423432)
hom_test_polysph(X = X, d = d, labels = hippocampus$ids_labs, type = "jsd",
                 h = 4 * h_rot, kernel = 3, kernel_type = 2, k = 100,
                 B = 100, show_prog = FALSE)

## ----citation, eval = FALSE---------------------------------------------------
# citation("polykde")