Nothing
exp_map <- function(psi, v, wnorm = l2_norm){
v_norm <- wnorm(v)
expgam <- cos(v_norm) * psi + sin(v_norm) * v / v_norm
return(expgam)
}
l2_curvenorm <- function(psi, time=seq(0,1,length.out=ncol(psi))){
sqrt(trapz(time,apply(psi^2,2,sum)))
}
#' map square root of warping function to tangent space
#'
#'
#' @param Psi vector describing psi function at center of tangent space
#' @param psi vector describing psi function to map to tangent space
#'
#' @return A numeric array of the same length as the input array `psi` storing the
#' shooting vector of `psi`
#'
#' @keywords srvf alignment
#' @export
inv_exp_map<-function(Psi, psi){
ip <- inner_product(Psi, psi)
if(ip < -1){
ip = -1
}else if(ip > 1){
ip = 1
}
theta <- acos(ip)
if (theta < 1e-10){
exp_inv = rep(0,length(psi))
} else {
exp_inv = theta / sin(theta) * (psi-cos(theta)*Psi)
}
return(exp_inv)
}
l2_norm<-function(psi, time=seq(0,1,length.out=length(psi))){
l2norm <- sqrt(trapz(time,psi*psi))
return(l2norm)
}
inner_product<-function(psi1, psi2, time=seq(0,1,length.out=length(psi1))){
ip <- trapz(time,psi1*psi2)
return(ip)
}
warp_q_gamma <- function(time, q, gam){
M = length(gam)
gam_dev = gradient(gam, 1/(M-1))
q_tmp = stats::approx(time,q,xout=(time[length(time)]-time[1])*gam +
time[1])$y*sqrt(gam_dev)
return(q_tmp)
}
randomGamma <- function(gam,num){
out = SqrtMean(gam)
mu = out$mu
psi = out$psi
vec = out$vec
K = stats::cov(t(vec))
out = svd(K)
s = out$d
U = out$u
n = 5
TT = nrow(vec)
vm = rowMeans(vec)
time <- seq(0,1,length.out=TT)
rgam = matrix(0,num,TT)
for (k in 1:num){
a = stats::rnorm(n)
v = rep(0,length(vm))
for (i in 1:n){
v = v + a[i]*sqrt(s[i])*U[,i]
}
psi <- exp_map(mu,v)
gam0 <- cumtrapz(time,psi*psi)
rgam[k,] = (gam0 - min(gam0))/(max(gam0)-min(gam0))
}
return(rgam)
}
#' SRVF transform of warping functions
#'
#' This function calculates the srvf of warping functions with corresponding
#' shooting vectors and finds the inverse of mean
#'
#' @param gam matrix (\eqn{N} x \eqn{M}) of \eqn{M} warping functions with \eqn{N} samples
#' @return `gamI` inverse of Karcher mean warping function
#'
#' @keywords srvf alignment
#' @references Srivastava, A., Wu, W., Kurtek, S., Klassen, E., Marron, J. S.,
#' May 2011. Registration of functional data using fisher-rao metric,
#' arXiv:1103.3817v2.
#' @references Tucker, J. D., Wu, W., Srivastava, A.,
#' Generative Models for Function Data using Phase and Amplitude Separation,
#' Computational Statistics and Data Analysis (2012), 10.1016/j.csda.2012.12.001.
#' @export
#' @examples
#' gamI <- SqrtMeanInverse(simu_warp$warping_functions)
SqrtMeanInverse <- function(gam){
TT = nrow(gam)
n = ncol(gam)
eps = .Machine$double.eps
time <- seq(0,1,length.out=TT)
psi = matrix(0,TT,n)
binsize <- mean(diff(time))
for (i in 1:n){
psi[,i] = sqrt(gradient(gam[,i],binsize))
}
# Find Direction
mu = rowMeans(psi)
stp <- .3
maxiter = 501
vec = matrix(0,TT,n)
lvm = rep(0,maxiter)
iter <- 1
for (i in 1:n){
vec[,i] <- inv_exp_map(mu, psi[,i])
}
vbar <- rowMeans(vec)
lvm[iter] <- l2_norm(vbar)
while (lvm[iter]>0.00000001 & iter<maxiter){
mu <- exp_map(mu, stp*vbar)
iter <- iter + 1
for (i in 1:n){
vec[,i] <- inv_exp_map(mu, psi[,i])
}
vbar <- rowMeans(vec)
lvm[iter] <- l2_norm(vbar)
}
gam_mu = cumtrapz(time, mu*mu)
gam_mu = (gam_mu - min(gam_mu))/(max(gam_mu)-min(gam_mu))
gamI = invertGamma(gam_mu)
return(gamI)
}
findkarcherinv <- function(warps, times, round = F){
m <- dim(warps)[1]
n <- dim(warps)[2]
psi.m <- matrix(0,m-1,n)
for(j in 1:n){psi.m[,j]<- sqrt(diff(warps[,j])/times)}
w <- apply(psi.m,1,mean)
mupsi <- w/sqrt(sum(w^2/(m-1)))
v.m <- matrix(0,m-1,n)
check <- 1
while(check > 0.01){
for (i in 1:n){
theta <- acos(sum(mupsi*psi.m[,i]/(m-1)))
v.m[,i] <- theta/sin(theta)*(psi.m[,i]-cos(theta)*mupsi)
}
vbar <- apply(v.m,1,mean)
check <- Enorm(vbar)/sqrt(m-1)
if (check>0){
mupsi.update <- cos(0.01*Enorm(vbar)/sqrt(m-1))*mupsi+sin(0.01*Enorm(vbar)/sqrt(m-1))*vbar/(Enorm(vbar)/sqrt(m-1))
}
else { mupsi.update <- cos(0.01*Enorm(vbar)/sqrt(m-1))*mupsi}
}
karcher.s <- 1+c(0,cumsum(mupsi.update^2)*times)
if(round){
invidy <- c(round(stats::approx(karcher.s,seq(1,(m-1)*times+1,times),method="linear",xout=1:((m-1)*times))$y),(m-1)*times+1)
}
else{
invidy <- c((stats::approx(karcher.s,seq(1,(m-1)*times+1,times),method="linear",xout=1:((m-1)*times))$y),(m-1)*times+1)
}
revscalevec <- sqrt(diff(invidy))
return(list(invidy = invidy,revscalevec = revscalevec))
}
#' map warping function to tangent space at identity
#'
#'
#' @param gam Either a numeric vector of a numeric matrix or a numeric array
#' specifying the warping functions
#' @param smooth Apply smoothing before gradient
#'
#' @return A numeric array of the same shape as the input array `gamma` storing the
#' shooting vectors of `gamma` obtained via finite differences.
#'
#' @keywords srvf alignment
#' @export
gam_to_v<-function(gam, smooth=TRUE){
if (ndims(gam) == 0){
TT = length(gam)
eps = .Machine$double.eps
time <- seq(0,1,length.out=TT)
binsize <- mean(diff(time))
psi = rep(0,TT)
if (smooth) {
tmp.spline <- stats::smooth.spline(gam)
g <- stats::predict(tmp.spline, deriv = 1)$y / binsize
g[g<0] = 0
psi = sqrt(g)
} else {
psi = sqrt(gradient(gam,binsize))
}
mu = rep(1,TT)
vec <- inv_exp_map(mu, psi)
} else {
TT = nrow(gam)
n = ncol(gam)
eps = .Machine$double.eps
time <- seq(0,1,length.out=TT)
binsize <- mean(diff(time))
psi = matrix(0,TT,n)
if (smooth) {
g <- matrix(0, TT, n)
for (i in 1:n) {
tmp.spline <- stats::smooth.spline(gam[,i])
g[, i] <- stats::predict(tmp.spline, deriv = 1)$y / binsize
g[g[,i]<0, i] = 0
psi[,i] = sqrt(g[, i])
}
} else {
for (i in 1:n){
psi[,i] = sqrt(gradient(gam[,i],binsize))
}
}
mu = rep(1,TT)
vec = matrix(0,TT,n)
for (i in 1:n){
vec[,i] <- inv_exp_map(mu, psi[,i])
}
}
return(vec)
}
#' map shooting vector to warping function at identity
#'
#'
#' @param v Either a numeric vector of a numeric matrix or a numeric array
#' specifying the shooting vectors
#'
#' @return A numeric array of the same shape as the input array `v` storing the
#' shooting vectors of `v` obtained via finite differences.
#'
#' @keywords srvf alignment
#' @export
v_to_gam<-function(v){
if (ndims(v) == 0){
TT = length(v)
time <- seq(0,1,length.out=TT)
mu = rep(1,TT)
psi <- exp_map(mu,v)
gam0 <- cumtrapz(time,psi*psi)
gam <- (gam0 - min(gam0))/(max(gam0)-min(gam0))
} else {
TT = nrow(v)
n = ncol(v)
time <- seq(0,1,length.out=TT)
mu = rep(1,TT)
gam = matrix(0,TT,n)
for (i in 1:n){
psi <- exp_map(mu,v[,i])
gam0 <- cumtrapz(time,psi*psi)
gam[,i] <- (gam0 - min(gam0))/(max(gam0)-min(gam0))
}
}
return(gam)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.