Merge branch 'master' of github.com:jdtuck/fdasrvf_R

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: fdasrvf
 Type: Package
 Title: Elastic Functional Data Analysis
-Version: 2.0.2
-Date: 2023-05-12
+Version: 2.1.1
+Date: 2023-12-15
 Authors@R: c(
     person(given = "J. Derek", 
            family = "Tucker", 

NAMESPACE

@@ -13,6 +13,7 @@ S3method(predict,mlpcr)
 S3method(predict,pcr)
 S3method(summary,fdawarp)
 export(SqrtMean)
+export(SqrtMeanInverse)
 export(SqrtMedian)
 export(align_fPCA)
 export(bootTB)
@@ -29,11 +30,16 @@ export(elastic.distance)
 export(elastic.lpcr.regression)
 export(elastic.mlpcr.regression)
 export(elastic.pcr.regression)
+export(elastic_amp_change_ff)
+export(elastic_change_fpca)
+export(elastic_ph_change_ff)
 export(f_to_srvf)
 export(function_group_warp_bayes)
+export(gam_to_v)
 export(gauss_model)
 export(gradient)
 export(horizFPCA)
+export(inv_exp_map)
 export(invertGamma)
 export(jointFPCA)
 export(joint_gauss_model)
@@ -55,6 +61,7 @@ export(sample_shapes)
 export(smooth.data)
 export(srvf_to_f)
 export(time_warping)
+export(v_to_gam)
 export(vertFPCA)
 export(warp_f_gamma)
 export(warp_q_gamma)

NEWS.md

@@ -1,4 +1,16 @@
-# fdasrvf 2.0.2.900
+# fdasrvf 2.1.1
+* bugfixes
+
+# fdasrvf 2.1.0
+* added elastic change point functions
+* exposed `SqrtMeanInverse` and `inv_exp_map` to global
+* bugfixes
+
+# fdasrvf 2.0.3
+
+* exposed lam to curve functions
+* added gamma to shooting vector conversion functions
+* bugfixes
 
 # fdasrvf 2.0.2
 

R/BBridge.R

@@ -0,0 +1,9 @@
+BBridge <- function(x=0, y=0, t0=0, T=1, N=100){
+  if(T<= t0) stop("wrong times")
+  dt <- (T-t0)/N
+  t <- seq(t0, T, length=N+1)
+  X <- c(0,cumsum(stats::rnorm(N)*sqrt(dt)))
+  BB <- x + X - (t-t0)/(T-t0)*(X[N+1]-y+x)
+  X <- stats::ts(BB, start=t0,deltat=dt)
+  return(invisible(X))
+}

R/LongRunCovMatrix.R

@@ -0,0 +1,65 @@
+#' Long Run Covariance Matrix Estimation for Multivariate Time Series
+#'
+#' This function estimates the long run covariance matrix of a given multivariate data sample.
+#'
+#' @param mdobj A multivariate data object
+#' @param h The bandwidth parameter. It is strictly non-zero. Choosing the bandwidth parameter to be zero is identical
+#' to estimating covariance matrix assuming iid data.
+#' @param kern_type Kernel function to be used for the estimation of the long run covariance
+#' matrix. The choices are \code{c("BT", "PR", "SP", "FT")} which are respectively, bartlett, parzen, simple and flat-top kernels.
+#' By default the function uses a \code{"barlett"} kernel.
+#' @return Returns long run covariance matrix
+
+# this is for the computation of Long Run Variance of \Theta
+LongRunCovMatrix <- function(mdobj, h=0, kern_type = "bartlett"){
+  N = ncol(mdobj)
+  D = nrow(mdobj)
+  Kernel <- function(i, h) {
+    x = i/h
+    if (kern_type == "flat") {
+      return(1)
+    }
+    if (kern_type == "simple") {
+      return(0)
+    }
+    if (kern_type == "bartlett") {
+      return(1 - x)
+    }
+    if (kern_type == "flat_top") {
+      if (x < 0.1) {
+        return(1)
+      } else {
+        if (x >= 0.1 & x < 1.1) {
+          return(1.1 - x)
+        } else {
+          return(0)
+        }
+      }
+    }
+    if (kern_type == "parzen") {
+      if (x < 1/2) {
+        return(1 - 6 * x^2 + 6 * abs(x)^3)
+      } else {
+        return(2 * (1 - abs(x))^3)
+      }
+    }
+  }
+  D_mat = matrix(0, D, D)
+  cdata = mdobj
+  # Long Run Cov Est
+  for (k in 1:D) {
+    for (r in k:D) {
+      s = cdata[k, 1:N] %*% cdata[r, 1:N]
+      if (h > 0) {
+        for (i in 1:h) {
+          a = cdata[k, 1:(N - i)] %*% cdata[r,(i + 1):N]
+          a = a + cdata[r, 1:(N - i)] %*% cdata[k,(i + 1):N]
+          s = s + Kernel(i, h) * a
+        }
+      }
+      D_mat[k, r] = s
+      D_mat[r, k] = D_mat[k, r]
+    }
+  }
+  D_mat/N
+}

R/curve_functions.R

@@ -2,16 +2,16 @@
 calculatecentroid <- function(beta,returnlength = F){
   n = nrow(beta)
   T1 = ncol(beta)
-  
+
   betadot = apply(beta,1,gradient,1.0/(T1-1))
   betadot = t(betadot)
-  
+
   normbetadot = apply(betadot,2,pvecnorm,2)
   integrand = matrix(0, n, T1)
   for (i in 1:T1){
     integrand[,i] = beta[,i] * normbetadot[i]
   }
-  
+
   scale = trapz(seq(0,1,length.out=T1), normbetadot)
   centroid = apply(integrand,1,trapz,x = seq(0,1,length.out=T1))/scale
   if(returnlength)  return(list("length" = scale,"centroid" = centroid))
@@ -140,10 +140,12 @@ shift_f <- function(f, tau){
     n = nrow(f)
     T1 = ncol(f)
     fn = matrix(0, n, T1)
-    if (tau > 0){
+    if (tau == T1){
+        fn[,(T1-tau+1):T1] = f[,1:tau]
+    } else if (tau > 0){
         fn[,1:(T1-tau)] = f[,(tau+1):T1]
         fn[,(T1-tau+1):T1] = f[,1:tau]
-    } else if (tau ==0) {
+    } else if (tau == 0) {
       fn = f
     } else {
         t = abs(tau)+1
@@ -169,7 +171,11 @@ find_rotation_seed_coord <- function(beta1, beta2, mode="O"){
 
     for (ctr in 0:end_idx){
         if (mode=="C"){
-            beta2n = shift_f(beta2, scl*ctr)
+            if ((scl*ctr) <= end_idx){
+              beta2n = shift_f(beta2, scl*ctr)
+            } else {
+              break
+            }
         } else {
             beta2n = beta2
         }
@@ -220,7 +226,7 @@ find_rotation_seed_coord <- function(beta1, beta2, mode="O"){
 }
 
 
-find_rotation_seed_unqiue <- function(q1, q2, mode="O"){
+find_rotation_seed_unqiue <- function(q1, q2, mode="O", lam=0.0){
     n1 = nrow(q1)
     T1 = ncol(q1)
     scl = 4
@@ -247,7 +253,7 @@ find_rotation_seed_unqiue <- function(q1, q2, mode="O"){
             dim(q1i) = c(T1*n1)
             q2i = q2n
             dim(q2i) = c(T1*n1)
-            gam0 = .Call('DPQ', PACKAGE = 'fdasrvf', q1i, q2i, n1, T1, 0, 1, 0, rep(0,T1))
+            gam0 = .Call('DPQ', PACKAGE = 'fdasrvf', q1i, q2i, n1, T1, lam, 1, 0, rep(0,T1))
             gamI = invertGamma(gam0)
             gam = (gamI-gamI[1])/(gamI[length(gamI)]-gamI[1])
             beta2n = q_to_curve(q2n)

R/curve_karcher_mean.R

@@ -6,6 +6,7 @@
 #' @param mode Open ("O") or Closed ("C") curves
 #' @param rotated Optimize over rotation (default = T)
 #' @param scale Include scale (default = F)
+#' @param lambda A numeric value specifying the elasticity. Defaults to `0.0`.
 #' @param maxit maximum number of iterations
 #' @param ms string defining whether the Karcher mean ("mean") or Karcher median ("median") is returned (default = "mean")
 #' @return Returns a list containing \item{mu}{mean srvf}
@@ -28,7 +29,8 @@
 #' @examples
 #' out <- curve_karcher_mean(beta[, , 1, 1:2], maxit = 2)
 #' # note: use more shapes, small for speed
-curve_karcher_mean <- function (beta, mode = "O", rotated = T, scale = F, maxit = 20, ms = "mean")
+curve_karcher_mean <- function (beta, mode = "O", rotated = T, scale = F,
+                                lambda = 0.0, maxit = 20, ms = "mean")
 {
     if(ms!="mean"&ms!="median"){warning("ms must be either \"mean\" or \"median\". ms has been set to \"mean\"",immediate. = T)}
     if(ms!="median"){ms = "mean"}
@@ -75,7 +77,7 @@ curve_karcher_mean <- function (beta, mode = "O", rotated = T, scale = F, maxit
     cat("\nInitializing...\n")
     gam = matrix(0,T1,N)
     for (k in 1:N) {
-        out = find_rotation_seed_unqiue(mu,q[, , k],mode)
+        out = find_rotation_seed_unqiue(mu,q[, , k],mode,lambda)
         gam[,k] = out$gambest
     }
 
@@ -96,7 +98,7 @@ curve_karcher_mean <- function (beta, mode = "O", rotated = T, scale = F, maxit
         for (i in 1:N) {
             q1 = q[, , i]
 
-            out = find_rotation_seed_unqiue(mu,q1,mode)
+            out = find_rotation_seed_unqiue(mu,q1,mode,lambda)
             qn_t = out$q2best/sqrt(innerprod_q2(out$q2best,out$q2best))
 
             q1dotq2 = innerprod_q2(mu,qn_t)

R/curve_pair_align.R

@@ -4,16 +4,19 @@
 #'
 #' @param beta1 array describing curve 1 (n,T)
 #' @param beta2 array describing curve 2 (n,T)
+#' @param mode Open ("O") or Closed ("C") curves
 #' @return a list containing \item{beta2n}{aligned curve 2 to 1}
 #' \item{q2n}{aligned srvf 2 to 1}
 #' \item{gam}{warping function}
 #' \item{q1}{srvf of curve 1}
+#' \item{beta1}{centered curve 1}
+#' \item{beta2}{centered curve 2}
 #' @keywords srvf alignment
 #' @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.
 #' @export
 #' @examples
 #' out <- curve_pair_align(beta[, , 1, 1], beta[, , 1, 5])
-curve_pair_align <- function(beta1, beta2){
+curve_pair_align <- function(beta1, beta2, mode="O"){
     T1 = ncol(beta1)
     centroid1 = calculatecentroid(beta1)
     dim(centroid1) = c(length(centroid1),1)
@@ -23,13 +26,11 @@ curve_pair_align <- function(beta1, beta2){
     beta2 = beta2 - repmat(centroid2, 1, T1)
 
     q1 = curve_to_q(beta1)$q
+    out = find_rotation_seed_coord(beta1, beta2, mode)
+    gam = out$gambest
+    q2n = out$q2best
+    beta2n = out$Rbest %*% shift_f(beta2, out$tau)
+    beta2n = group_action_by_gamma_coord(beta2n, gam)
 
-    # optimize over SO(n) x Gamma using DP
-    out = reparam_curve(beta1, beta2)
-    beta2n = out$R %*% shift_f(beta2, out$tau)
-    gamI = invertGamma(out$gam)
-    beta2n = group_action_by_gamma_coord(beta2n, gamI)
-    q2n = curve_to_q(beta2n)$q
-
-    return(list(beta2n=out$beta2new, q2n=q2n, gam=gamI, q1=q1))
+    return(list(beta2n=beta2n, q2n=q2n, gam=gam, q1=q1, beta1=beta1, beta2=beta2))
 }

R/curve_srvf_align.R

@@ -6,6 +6,7 @@
 #' @param mode Open ("O") or Closed ("C") curves
 #' @param rotated Optimize over rotation (default = T)
 #' @param scale Include scale (default = F)
+#' @param lambda A numeric value specifying the elasticity. Defaults to `0.0`.
 #' @param maxit maximum number of iterations
 #' @param ms string defining whether the Karcher mean ("mean") or Karcher median ("median") is returned (default = "mean")
 #' @return Returns a list containing \item{betan}{aligned curves}
@@ -17,16 +18,18 @@
 #' @export
 #' @examples
 #' data("mpeg7")
-#' out = curve_srvf_align(beta[,,1,1:2],maxit=2) # note: use more shapes, small for speed
-curve_srvf_align <- function(beta, mode="O", rotated=T, scale = F, maxit=20, ms = "mean"){
+#' # note: use more shapes and iterations, small for speed
+#' out = curve_srvf_align(beta[,,1,1:2],maxit=2)
+curve_srvf_align <- function(beta, mode="O", rotated=T, scale = F, lambda = 0.0,
+                             maxit=20, ms = "mean"){
     if (mode=="C"){
       isclosed = TRUE
     }
     tmp = dim(beta)
     n = tmp[1]
     T1 = tmp[2]
     N = tmp[3]
-    out = curve_karcher_mean(beta, mode, rotated, scale, maxit, ms)
+    out = curve_karcher_mean(beta, mode, rotated, scale, lambda, maxit, ms)
     beta<-out$beta
     mu = out$mu
     betamean = out$betamean
@@ -35,24 +38,24 @@ curve_srvf_align <- function(beta, mode="O", rotated=T, scale = F, maxit=20, ms
 
     qn = array(0, c(n,T1,N))
     betan = array(0, c(n,T1,N))
-    rotmat = array(0, c(n,n,N)) 
-    gams = matrix(0, T1, N)     
- 
+    rotmat = array(0, c(n,n,N))
+    gams = matrix(0, T1, N)
+
     # align to mean
     for (ii in 1:N){
         q1 = q[,,ii]
         beta1 = beta[,,ii]
 
-        out = find_rotation_seed_unqiue(mu,q1,mode)
-        gams[,ii] = out$gambest 
+        out = find_rotation_seed_unqiue(mu,q1,mode,lambda)
+        gams[,ii] = out$gambest
         beta1 = out$Rbest%*%beta1
         beta1n = group_action_by_gamma_coord(beta1, out$gambest)
         q1n = curve_to_q(beta1n)$q
 
         out = find_best_rotation(mu, q1n)
         qn[,,ii] = out$q2new
         betan[,,ii] = out$R%*%beta1n
-        rotmat[,,ii] = out$R  
+        rotmat[,,ii] = out$R
     }
     return(list(betan=betan, qn=qn, betamean=betamean, q_mu=mu, rotmat = rotmat,gams = gams,v=v))
 }