Merge pull request #47 from astamm/karcher_mean

Karcher mean and draft proposal of refactoring

NAMESPACE

@@ -22,6 +22,7 @@ export(SqrtMedian)
 export(align_fPCA)
 export(bootTB)
 export(calc_shape_dist)
+export(curve2srvf)
 export(curve_boxplot)
 export(curve_depth)
 export(curve_dist)
@@ -32,6 +33,8 @@ export(curve_pair_align)
 export(curve_principal_directions)
 export(curve_srvf_align)
 export(curve_to_q)
+export(discrete2curve)
+export(discrete2warping)
 export(elastic.depth)
 export(elastic.distance)
 export(elastic.lpcr.regression)
@@ -44,6 +47,13 @@ export(f_to_srvf)
 export(function_group_warp_bayes)
 export(gam_to_v)
 export(gauss_model)
+export(get_distance_matrix)
+export(get_hilbert_sphere_distance)
+export(get_identity_warping)
+export(get_l2_distance)
+export(get_l2_inner_product)
+export(get_shape_distance)
+export(get_warping_distance)
 export(gradient)
 export(horizFPCA)
 export(inv_exp_map)
@@ -69,12 +79,16 @@ export(rgam)
 export(sample_shapes)
 export(shape_CI)
 export(smooth.data)
+export(srvf2curve)
 export(srvf_to_f)
 export(time_warping)
+export(to_hilbert_sphere)
 export(v_to_gam)
 export(vertFPCA)
+export(warp_curve)
 export(warp_f_gamma)
 export(warp_q_gamma)
+export(warp_srvf)
 importFrom(Rcpp,sourceCpp)
 importFrom(coda,effectiveSize)
 importFrom(foreach,"%dopar%")

R/calc_shape_dist.R

@@ -1,9 +1,13 @@
 #' Elastic Shape Distance
 #'
 #' Calculates the elastic shape distance between two curves `beta1` and `beta2`.
-#' If the input curves are describing multidimensional functional data, then
-#' most of the time the user should set `rotation == FALSE`, `scale == FALSE`
-#' and `mode == "O"`.
+#'
+#' Distances are computed between the SRVFs of the original curves. Hence, they
+#' are intrinsically invariant to position. The user can then choose to make
+#' distances invariant to rotation and scaling by setting `rotation` and `scale`
+#' accordingly. Distances can also be made invariant to reparametrization by
+#' setting `alignment = TRUE`, in which case curves are aligned using an
+#' appropriate action of the diffeomorphism group on the space of SRVFs.
 #'
 #' @param beta1 A numeric matrix of shape \eqn{L \times M} specifying an
 #'   \eqn{L}-dimensional curve evaluated on \eqn{M} sample points.
@@ -25,6 +29,9 @@
 #'   about the actual length of the SRVFs in the metric when using normalized
 #'   SRVFs to achieve scale invariance. This only applies if `scale == TRUE`.
 #'   Defaults to `FALSE`.
+#' @param lambda A numeric value specifying the weight of a penalty term that
+#'   constraints the warping function to be not too far from the identity.
+#'   Defaults to `0.0`.
 #'
 #' @return A list with the following components:
 #' - `d`: the amplitude (geodesic) distance;
@@ -49,6 +56,8 @@
 #' @references Kurtek, S., Srivastava, A., Klassen, E., and Ding, Z. (2012),
 #'   “Statistical Modeling of Curves Using Shapes and Related Features,” Journal
 #'   of the American Statistical Association, 107, 1152–1165.
+#' @references Srivastava, A., Klassen, E. P. (2016). Functional and shape
+#'   data analysis, 1. New York: Springer.
 #'
 #' @export
 #' @examples
@@ -58,7 +67,12 @@ calc_shape_dist <- function(beta1, beta2,
                             alignment = TRUE,
                             rotation = TRUE,
                             scale = TRUE,
-                            include.length = FALSE) {
+                            include.length = FALSE,
+                            lambda = 0.0) {
+  if (mode == "C" && !scale)
+    cli::cli_abort("Closed curves are currently handled only on the L2
+                   hypersphere. Please set `scale = TRUE`.")
+
   srvf1 <- curve_to_srvf(beta1, scale = scale)
   srvf2 <- curve_to_srvf(beta2, scale = scale)
 
@@ -67,7 +81,9 @@ calc_shape_dist <- function(beta1, beta2,
     mode = mode,
     alignment = alignment,
     rotation = rotation,
-    scale = scale
+    scale = scale,
+    lambda = lambda,
+    include_length = include.length
   )
 
   list(

R/curve_dist.R

@@ -3,6 +3,8 @@
 #' Computes the pairwise distance matrix between a set of curves using the
 #' elastic shape distance as computed by [`calc_shape_dist()`].
 #'
+#' @inherit calc_shape_dist details
+#'
 #' @param beta A numeric array of shape \eqn{L \times M \times N} specifying the
 #'   set of \eqn{N} curves of length \eqn{M} in \eqn{L}-dimensional space.
 #' @inheritParams calc_shape_dist
@@ -15,6 +17,17 @@
 #'   the amplitude and phase distances, respectively.
 #' @export
 #'
+#' @keywords distances
+#'
+#' @references Srivastava, A., Klassen, E., Joshi, S., Jermyn, I., (2011). Shape
+#'   analysis of elastic curves in euclidean spaces. Pattern Analysis and
+#'   Machine Intelligence, IEEE Transactions on 33 (7), 1415-1428.
+#' @references Kurtek, S., Srivastava, A., Klassen, E., and Ding, Z. (2012),
+#'   “Statistical Modeling of Curves Using Shapes and Related Features,” Journal
+#'   of the American Statistical Association, 107, 1152–1165.
+#' @references Srivastava, A., Klassen, E. P. (2016). Functional and shape
+#'   data analysis, 1. New York: Springer.
+#'
 #' @examples
 #' out <- curve_dist(beta[, , 1, 1:4])
 curve_dist <- function(beta,
@@ -23,7 +36,12 @@ curve_dist <- function(beta,
                        rotation = FALSE,
                        scale = FALSE,
                        include.length = FALSE,
+                       lambda = 0.0,
                        ncores = 1L) {
+  if (mode == "C" && !scale)
+    cli::cli_abort("Closed curves are currently handled only on the L2
+                   hypersphere. Please set `scale = TRUE`.")
+
   navail <- max(parallel::detectCores() - 1, 1)
 
   if (ncores > navail) {
@@ -46,8 +64,9 @@ curve_dist <- function(beta,
   M <- dims[2]
   N <- dims[3]
 
+  srvfs <- lapply(1:N, \(n) curve_to_srvf(beta[, , n], scale = scale))
   for (n in 1:N)
-    beta[ , , n] <- curve_to_srvf(beta[, , n], scale = scale)$q
+    beta[ , , n] <- srvfs[[n]]$q
 
   K <- N * (N - 1) / 2
 
@@ -61,18 +80,41 @@ curve_dist <- function(beta,
     q1 <- beta[, , i + 1]
     q2 <- beta[, , j + 1]
 
-    out <- find_rotation_seed_unique(
+    norm_ratio <- 1
+    if (scale)
+      norm_ratio <- srvfs[[i + 1]]$qnorm / srvfs[[j + 1]]$qnorm
+
+    res12 <- find_rotation_seed_unique(
       q1, q2,
       mode = mode,
       alignment = alignment,
       rotation = rotation,
-      scale = scale
+      scale = scale,
+      norm_ratio = norm_ratio,
+      lambda = lambda
+    )
+
+    res21 <- find_rotation_seed_unique(
+      q2, q1,
+      mode = mode,
+      alignment = alignment,
+      rotation = rotation,
+      scale = scale,
+      norm_ratio = norm_ratio,
+      lambda = lambda
     )
+
+    if (res12$d < res21$d)
+      res <- res12
+    else
+      res <- res21
+
     if (alignment)
-      dx <- phase_distance(out$gambest)
+      dx <- phase_distance(res$gambest)
     else
       dx <- 0
-    matrix(c(out$d, dx), ncol = 1)
+
+    matrix(c(res$d, dx), ncol = 1)
   }
 
   Da <- out[1, ]

R/curve_functions.R

@@ -246,9 +246,18 @@ find_rotation_seed_unique <- function(q1, q2,
                                       alignment = TRUE,
                                       rotation = TRUE,
                                       scale = TRUE,
+                                      norm_ratio = 1.0,
                                       lambda = 0.0) {
   L <- nrow(q1)
   M <- ncol(q1)
+
+  # Variables for DPQ2 algorithm
+  grd <- seq(0, 1, length.out = M)
+  nbhd_dim <- 7L
+  Gvec <- rep(0, M)
+  Tvec <- rep(0, M)
+  size <- 0
+
   scl <- 4
   minE <- Inf
 
@@ -280,9 +289,16 @@ find_rotation_seed_unique <- function(q1, q2,
       }
       dim(q1i) <- M * L
       dim(q2ni) <- M * L
-      gam0 <- .Call('DPQ', PACKAGE = 'fdasrvf', q1i, q2ni,
-                    L, M, lambda, 1, 0, rep(0, M))
-      gamI <- invertGamma(gam0)
+
+      ret <- .Call(
+        "DPQ2", PACKAGE = "fdasrvf",
+        q1i, grd, q2ni, grd, L, M, M, grd, grd, M,
+        M, Gvec, Tvec, size, lambda, nbhd_dim
+      )
+      Gvec <- ret$G[1:ret$size]
+      Tvec <- ret$T[1:ret$size]
+      gamI <- stats::approx(Tvec, Gvec, xout = grd)$y
+
       gam <- (gamI - gamI[1]) / (gamI[length(gamI)] - gamI[1])
       q2new <- group_action_by_gamma(q2n, gam, scale = scale)
 
@@ -293,7 +309,7 @@ find_rotation_seed_unique <- function(q1, q2,
       gam <- seq(0, 1, length.out = M)
     }
 
-    Ec <- amplitude_distance(q1, q2new, scale = scale)
+    Ec <- amplitude_distance(q1, q2new, scale = scale, norm_ratio = norm_ratio)
 
     if (Ec < minE) {
       Rbest <- Rbest
@@ -336,9 +352,9 @@ find_rotation_and_seed_q <- function(q1,q2){
 group_action_by_gamma <- function(q, gamma, scale = TRUE) {
   L <- nrow(q)
   M <- ncol(q)
+  grd <- seq(0, 1, length.out = M)
   gammadot <- gradient(gamma, 1.0 / M)
   qn <- matrix(nrow = L, ncol = M)
-  grd <- seq(0, 1, length.out = M)
 
   for (l in 1:L)
     qn[l, ] <- stats::spline(grd, q[l, ], xout = gamma)$y * sqrt(gammadot)
@@ -721,13 +737,20 @@ match_f2_to_f1 <- function(srvf1, srvf2, beta2,
                            mode = "O",
                            alignment = TRUE,
                            rotation = FALSE,
-                           scale = FALSE) {
+                           scale = FALSE,
+                           include_length = FALSE,
+                           lambda = 0.0) {
+  norm_ratio <- 1
+  if (include_length)
+    norm_ratio <- srvf1$qnorm / srvf2$qnorm
   out <- find_rotation_seed_unique(
     srvf1$q, srvf2$q,
     mode = mode,
     alignment = alignment,
     rotation = rotation,
-    scale = scale
+    scale = scale,
+    lambda = lambda,
+    norm_ratio = norm_ratio
   )
 
   # Compute amplitude distance
@@ -743,7 +766,7 @@ match_f2_to_f1 <- function(srvf1, srvf2, beta2,
     qscale <- srvf1$qnorm / srvf2$qnorm
   betascale <- qscale^2
 
-  if (mode == "C") {
+  if (mode == "C" && scale) {
     beta2n <- q_to_curve(out$q2best, scale = qscale)
   } else {
     beta2n <- beta2 - srvf2$centroid

R/fdasrvf-package.R

@@ -61,4 +61,3 @@
 #' @importFrom foreach %dopar%
 #' @aliases fdasrvf fdasrvf-package
 "_PACKAGE"
-NULL

R/invertGamma.R

@@ -14,10 +14,10 @@
 #' @export
 #' @examples
 #' out <- invertGamma(simu_warp$warping_functions[, 1])
-invertGamma <- function(gam){
-    N = length(gam)
-    x = (0:(N-1))/(N-1)
-    gamI = stats::approx(gam,x,xout=x)$y
+invertGamma <- function(gam) {
+    N <- length(gam)
+    x <- seq(0, 1, len = N)
+    gamI = stats::approx(gam, x, xout = x)$y
     gamI[N] = 1
     gamI = gamI/gamI[N]
     return(gamI)