# R/shapes.R In shapes: Statistical Shape Analysis

```#-----------------------------------------------------------------------
#
# Statistical shape analysis routines
# written by Ian Dryden in R  (see http://cran.r-project.org)
# (c) Ian Dryden
#     UoN, FIU. version 1.2.7
#                    2003-2023
#
# Includes contributions by many other authors, including
#   Mohammad Faghihi, Kwang-Rae Kim, Alfred Kume,
#   Gregorio Quintana-Orti, Amelia Simo.
#
###########################################################################

tangentcoords.partial.inv = function(v, p, R)
{
return(matrix(sqrt(1 - sum(v^2)) * c(p) + v, nrow = nrow(p)) %*% t(R))
}

preshape2shape = function(z)
{
k = nrow(z) + 1
H = defh(k - 1)
return(t(H) %*% z)
}

plot3darcs<-function(x,pcno=1,c=1,nn=100,boundary.data=TRUE,view.theta=0,view.phi=0,type="pnss"){
# points along principal arcs

pns.out <- x
k <- pns.out\$GPAout\$k
m <- pns.out\$GPAout\$m
n.pc <- dim(pns.out\$resmat)[1]

npts = 100
arc1 = t(get.prinarc(resmat = pns.out\$resmat, PNS = pns.out\$PNS,
arc = 1, n = npts, boundary.data = boundary.data))
arc2 = t(get.prinarc(resmat = pns.out\$resmat, PNS = pns.out\$PNS,
arc = 2, n = npts, boundary.data = boundary.data))
arc3 = t(get.prinarc(resmat = pns.out\$resmat, PNS = pns.out\$PNS,
arc = 3, n = npts, boundary.data = boundary.data))
PNSmean = pns.out\$PNS\$mean
GPAout = pns.out\$GPAout
{
#     cat("stdev of PNS1 score:", round(sd(pns.out\$resmat[1,
#         ]), 4), "\n")
#    cat("stdev of PNS2 score:", round(sd(pns.out\$resmat[2,
#        ]), 4), "\n")
#    cat("stdev of PNS3 score:", round(sd(pns.out\$resmat[3,
#        ]), 4), "\n")
}
rng = c * sd(pns.out\$resmat[1, ])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)], 0, seq(0,
rng, length = nn + 1)[-1])
lu.arc1 = t(get.prinarc.value(PNS = pns.out\$PNS, arc = 1,
res = val))
rng = c * sd(pns.out\$resmat[2, ])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)], 0, seq(0,
rng, length = nn + 1)[-1])
lu.arc2 = t(get.prinarc.value(PNS = pns.out\$PNS, arc = 2,
res = val))
rng = c * sd(pns.out\$resmat[3, ])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)], 0, seq(0,
rng, length = nn + 1)[-1])
lu.arc3 = t(get.prinarc.value(PNS = pns.out\$PNS, arc = 3,
res = val))
scores.arc1 = sphere2pcscore(x = arc1)
scores.arc2 = sphere2pcscore(x = arc2)
scores.arc3 = sphere2pcscore(x = arc3)
scores.PNSmean = sphere2pcscore(x = t(PNSmean))
scores.lu.arc1 = sphere2pcscore(x = lu.arc1)
scores.lu.arc2 = sphere2pcscore(x = lu.arc2)
scores.lu.arc3 = sphere2pcscore(x = lu.arc3)
U1 = matrix(0, npts, nrow(GPAout\$pcar))
U2 = matrix(0, npts, nrow(GPAout\$pcar))
U3 = matrix(0, npts, nrow(GPAout\$pcar))
for (i in 1:npts) {
for (j in 1:n.pc) {
U1[i, ] = U1[i, ] + scores.arc1[i, j] * GPAout\$pcar[,
j]
U2[i, ] = U2[i, ] + scores.arc2[i, j] * GPAout\$pcar[,
j]
U3[i, ] = U3[i, ] + scores.arc3[i, j] * GPAout\$pcar[,
j]
}
}
U.mean = matrix(0, 1, nrow(GPAout\$pcar))
for (j in 1:n.pc) {
U.mean = U.mean + scores.PNSmean[j] * GPAout\$pcar[, j]
}
tan.lu.arc1 = matrix(0, nrow(lu.arc1), nrow(GPAout\$pcar))
tan.lu.arc2 = matrix(0, nrow(lu.arc2), nrow(GPAout\$pcar))
tan.lu.arc3 = matrix(0, nrow(lu.arc3), nrow(GPAout\$pcar))
for (i in 1:nrow(lu.arc1)) {
for (j in 1:n.pc) {
tan.lu.arc1[i, ] = tan.lu.arc1[i, ] + scores.lu.arc1[i,
j] * GPAout\$pcar[, j]
tan.lu.arc2[i, ] = tan.lu.arc2[i, ] + scores.lu.arc2[i,
j] * GPAout\$pcar[, j]
tan.lu.arc3[i, ] = tan.lu.arc3[i, ] + scores.lu.arc3[i,
j] * GPAout\$pcar[, j]
}
}
shapes.arc1 = array(NA, c(k, m, npts))
shapes.arc2 = array(NA, c(k, m, npts))
shapes.arc3 = array(NA, c(k, m, npts))
H = defh(k - 1)
for (i in 1:npts) {
# to convert from in expo map to partial tangent coords and then to icon configuration
rho<-Enorm(U1[i,])
shapes.arc1[, , i] = preshape2shape(tangentcoords.partial.inv(v = U1[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
rho<-Enorm(U2[i,])
shapes.arc2[, , i] = preshape2shape(tangentcoords.partial.inv(v = U2[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
rho<-Enorm(U3[i,])
shapes.arc3[, , i] = preshape2shape(tangentcoords.partial.inv(v = U3[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
}
rho<-Enorm(U.mean)
shapes.PNSmean = preshape2shape(tangentcoords.partial.inv(v = U.mean*sin(rho)/rho,
p = H %*% GPAout\$mshape, R = diag(m)))
shapes.lu.arc1 = array(NA, c(k, m, nrow(lu.arc1)))
shapes.lu.arc2 = array(NA, c(k, m, nrow(lu.arc2)))
shapes.lu.arc3 = array(NA, c(k, m, nrow(lu.arc3)))
for (i in 1:nrow(lu.arc1)) {
rho<-Enorm(tan.lu.arc1[i,])
shapes.lu.arc1[, , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc1[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
rho<-Enorm(tan.lu.arc2[i,])
shapes.lu.arc2[, , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc2[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
rho<-Enorm(tan.lu.arc3[i,])
shapes.lu.arc3[, , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc3[i,
]*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))
}
h <- defh(k - 1)
zero <- matrix(0, k - 1, k)
H <- cbind(h, zero, zero)
H1 <- cbind(zero, h, zero)
H2 <- cbind(zero, zero, h)
H <- rbind(H, H1, H2)
if (dim(GPAout\$pcar)[1] == (3 * (k - 1))) {
pcarot <- (t(H) %*% GPAout\$pcar)
GPAout\$pcar <- pcarot
}
if (pcno == 1) {
shapes.lu.arc <- shapes.lu.arc1
}
if (pcno == 2) {
shapes.lu.arc <- shapes.lu.arc2
}
if (pcno == 3) {
shapes.lu.arc <- shapes.lu.arc3
}
if (type == "pca") {
open3d()
par3d(windowRect = c(20, 30, 800, 800))
view3d(view.theta, view.phi)
lines3d(GPAout\$mshape, col = rainbow(k), lwd = 5)
pcu <- GPAout\$mshape + c * GPAout\$pcasd[pcno] * cbind(GPAout\$pcar[1:k,
pcno], GPAout\$pcar[(k + 1):(2 * k), pcno], GPAout\$pcar[(2 *
k + 1):(3 * k), pcno])
pcl <- GPAout\$mshape - c * GPAout\$pcasd[pcno] * cbind(GPAout\$pcar[1:k,
pcno], GPAout\$pcar[(k + 1):(2 * k), pcno], GPAout\$pcar[(2 *
k + 1):(3 * k), pcno])
for (j in 1:k) {
lines3d(rbind(pcl[j, ], pcu[j, ]), col = rainbow(k)[j])
if (j > 1) {
lines3d(rbind(pcu[j - 1, ], pcu[j, ]), col = "black")
lines3d(rbind(pcl[j - 1, ], pcl[j, ]), col = "grey")
}
}
}
if (type == "pnss") {
open3d()
par3d(windowRect = c(20, 30, 800, 800))
view3d(view.theta, view.phi)
plot3d(shapes.PNSmean, type = "s", col = rainbow(k),
lines3d(shapes.PNSmean, lwd = 5, col = rainbow(k))
for (i in 1:k) {
lines3d(t(shapes.lu.arc[i, , ]), col = rainbow(k)[i],
lwd = 1, lty = 2)
color = "black")
if (i > 1) {
lines3d((shapes.lu.arc[(i - 1):i, , 1]), col = "black")
}
color = "grey")
if (i > 1) {
lines3d((shapes.lu.arc[(i - 1):i, , 201]), col = "grey")
}
}
}
out <- list(PNSmean = 0, lu.arc = 0)
out\$PNSmean <- shapes.PNSmean
out\$lu.arc <- shapes.lu.arc
out
}

########
pnss3d<-
function (x, sphere.type = "seq.test", alpha = 0.1, R = 100,
nlast.small.sphere = 1, n.pc = "Full", output=TRUE)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n.pc =="Full" ) {
n.pc=m*k-m*(m-1)/2-m
}
if (m==2){
tem1 <- array( 0, c(k,3,n) )
tem1[,1:2,]<-x
x<-tem1
m<-3
}

#if (n < ((k - 1) * m)) {
#  print("Note: n < (k - 1) * m.")
#  jj<- round( (k-1)*m/n + 0.5)
#  print("Adding extra copies of the data")
#  tem<- array(0,c(k,m,jj*n))
#  tem[,,1:n]<-x
#  for (i in 2:jj){
#    for (j in 1:n){
#    tem[,,(i-1)*n+ j ]<-x[,,j] + 0*matrix( rnorm(k*m), k,m)
#    }
#  }
#  x<-tem
#}
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]

out = pc2sphere2(x = x, n.pc = n.pc, output=output)
spheredata = t(out\$spheredata)
GPAout = out\$GPAout
pns.out = pns(x = spheredata, sphere.type = sphere.type,
alpha = alpha, R = R, nlast.small.sphere = nlast.small.sphere, output=output)
pns.out\$percent = pns.out\$percent * sum(GPAout\$percent[1:n.pc])/100
if (output){
print("PNS percent explained")
cat(c(round(pns.out\$percent,2),"\n"))
print("PCA percent explained")
cat(c(round(GPAout\$percent,2),"\n"))
}
pns.out\$GPAout = GPAout
pns.out\$spheredata = spheredata
return(pns.out)
}

pc2sphere2<-function (x, n.pc, output=TRUE)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]

GPAout = procGPA(x = x, scale = TRUE, reflect = FALSE, tol1=1e-8,tangentcoords = "partial",
distances = TRUE)
if (output){
cat("First ", n.pc, " principal components explain ", round(sum(GPAout\$percent[1:n.pc]),2),
"% of total variance. \n", sep = "")
}
H = defh(k - 1)
X.hat = H %*% GPAout\$mshape
S = array(NA, c(k - 1, m, n))
for (i in 1:n) {
S[, , i] = H %*% GPAout\$rotated[, , i]
}
T.c = GPAout\$tan #- apply(GPAout\$tan, 1, mean)
out = pcscore2sphere2(n.pc = n.pc, X.hat = X.hat, S = S, Tan = T.c,
V = GPAout\$pcar)
return(list(spheredata = out, GPAout = GPAout))
}

backfit <- function( scores, x , type="pnss", size=1){
npc <- length(scores)

if (type=="pnss"){
PNS.object<-x
PNS<-PNS.object\$PNS
GPAout<-PNS.object\$GPAout
z1 <- PNSe2s(matrix(scores,npc,1),PNS)
pcscores<-c(sphere2pcscore(x=t(z1)))
#note the PC scores are from the inverse exponential map tangent coordinates
mu <- GPAout\$mshape
k<-dim(mu)[1]
m<-dim(mu)[2]
H = defh(k - 1)
U<- GPAout\$pcar[,1]*0
for (j in 1:npc) {
U = U + pcscores[j] * GPAout\$pcar[, j]
}
# to convert from in expo map to partial tangent coords and then to icon configuration
rho<-Enorm(U)
xout<-preshape2shape(tangentcoords.partial.inv(v = U*sin(rho)/rho, p = H %*% GPAout\$mshape, R = diag(m)))*size
}

if (type=="pca"){
GPAout<- x
pcscores<-scores #assume partial tangent coordinates
mu <- GPAout\$mshape
k<-dim(mu)[1]
m<-dim(mu)[2]
H = defh(k - 1)
U<- GPAout\$pcar[,1]*0
for (j in 1:npc) {
U = U + pcscores[j] * GPAout\$pcar[, j]
}
xout<-preshape2shape(tangentcoords.partial.inv(v = U, p = H %*% GPAout\$mshape, R = diag(m)))*size
}

xout
}

#==================================================================================
# PNS  The Principal Nested Spheres code (PNS) for spheres and shapes has
# been written by Kwang-Rae Kim, and builds closely on the original matlab
# code for PNS by Sungkyu Jung
#==================================================================================

#==================================================================================
pns = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 1, output=TRUE)
{
n = ncol(x)
k = nrow(x)
PNS = list()

if (abs(sum(apply(x ^ 2, 2, sum)) - n) > 1e-8)
{
stop("Error: Each column of x should be a unit vector, ||x[ , i]|| = 1.")
}
svd.x = svd(x, nu = nrow(x))
uu = svd.x\$u
maxd = which(svd.x\$d < 1e-15)[1]
if (is.na(maxd) | k > n)
{
maxd = min(k, n) + 1
}
nullspdim = k - maxd + 1
d = k - 1
if (output){
cat("Message from pns() : dataset is on ", d, "-sphere. \n", sep = "")
}
if (nullspdim > 0)
{
if (output){
cat(" .. found null space of dimension ",
nullspdim,
", to be trivially reduced. \n",
sep = "")
}
}

if (d==2){
PNS\$spherePNS<-t(x)
}

resmat = matrix(NA, d, n)
orthaxis = list()
orthaxis[[d - 1]] = NA
dist = rep(NA, d - 1)
pvalues = matrix(NA, d - 1, 2)
ratio = rep(NA, d - 1)
currentSphere = x
if (nullspdim > 0)
{
for (i in 1:nullspdim)
{
oaxis = uu[, ncol(uu) - i + 1]
r = pi / 2
pvalues[i,] = c(NaN, NaN)
res = acos(t(oaxis) %*% currentSphere) - r
orthaxis[[i]] = oaxis
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(oaxis) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
uu = rotMat(oaxis) %*% uu
uu = uu[1:(k - i),] / repmat(matrix(sqrt(1 - uu[nrow(uu),] ^ 2), nrow = 1), k - i, 1)
if (output){
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"NULL space \n",
sep = "")
}
}
}
if (sphere.type == "seq.test")
{
if (output){
cat(" .. sequential tests with significance level ",
alpha,
"\n",
sep = "")
}
isIsotropic = FALSE
for (i in (nullspdim + 1):(d - 1))
{
if (!isIsotropic)
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center.s = sp\$center
r.s = sp\$r
resSMALL = acos(t(center.s) %*% currentSphere) - r.s
sp = getSubSphere(x = currentSphere, geodesic = "great")
center.g = sp\$center
r.g = sp\$r
resGREAT = acos(t(center.g) %*% currentSphere) - r.g
pval1 = LRTpval(resGREAT, resSMALL, n)
pvalues[i, 1] = pval1
if (pval1 > alpha)
{
center = center.g
r = r.g
pvalues[i, 2] = NA
if (output){
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, p(LRT) = ",
pval1,
"\n",
sep = ""
)
}
} else {
pval2 = vMFtest(currentSphere, R)
pvalues[i, 2] = pval2
if (pval2 > alpha)
{
center = center.g
r = r.g
if (output){
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, p(LRT) = ",
pval1,
", p(vMF) = ",
pval2,
"\n",
sep = ""
)
}
isIsotropic = TRUE
} else {
center = center.s
r = r.s
if (output){
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by SMALL sphere, p(LRT) = ",
pval1,
", p(vMF) = ",
pval2,
"\n",
sep = ""
)
}
}
}
} else if (isIsotropic) {
sp = getSubSphere(x = currentSphere, geodesic = "great")
center = sp\$center
r = sp\$r
if (output){
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, restricted by testing vMF distn",
"\n",
sep = ""
)
}
pvalues[i, 1] = NA
pvalues[i, 2] = NA
}
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
###########
if (nrow(currentSphere) == 3)
{
PNS\$spherePNS = t(currentSphere)
}

if (nrow(currentSphere) == 2)
{
PNS\$circlePNS = t(cur.proj)
}
#############################
}
} else if (sphere.type == "BIC") {
if (output){
cat(" .. with BIC \n")
}
for (i in (nullspdim + 1):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center.s = sp\$center
r.s = sp\$r
resSMALL = acos(t(center.s) %*% currentSphere) - r.s
sp = getSubSphere(x = currentSphere, geodesic = "great")
center.g = sp\$center
r.g = sp\$r
resGREAT = acos(t(center.g) %*% currentSphere) - r.g
BICsmall = n * log(mean(resSMALL ^ 2)) + (d - i + 1 + 1) * log(n)
BICgreat = n * log(mean(resGREAT ^ 2)) + (d - i + 1) * log(n)
if (output){
cat("BICsm: ", BICsmall, ", BICgr: ", BICgreat, "\n", sep = "")
}
if (BICsmall > BICgreat)
{
center = center.g
r = r.g
if (output){
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"GREAT sphere, BIC \n",
sep = "")
}
} else {
center = center.s
r = r.s
if (output){
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"SMALL sphere, BIC \n",
sep = "")
}
}
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res

cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)

NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)

###########
if (nrow(currentSphere) == 3)
{
PNS\$spherePNS = t(currentSphere)
}

if (nrow(currentSphere) == 2)
{
PNS\$circlePNS = t(cur.proj)
}
#############################

}
} else if (sphere.type == "small" | sphere.type == "great") {
pvalues = NaN
for (i in (nullspdim + 1):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = sphere.type)
center = sp\$center
r = sp\$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res

cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)

NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)

###########
if (nrow(currentSphere) == 3)
{
PNS\$spherePNS = t(currentSphere)
}

if (nrow(currentSphere) == 2)
{
PNS\$circlePNS = t(cur.proj)
}
#############################

}
} else if (sphere.type == "bi.sphere") {
if (nlast.small.sphere < 0)
{
cat("!!! Error from pns(): \n")
cat("!!! nlast.small.sphere should be >= 0. \n")
return(NULL)
}
mx = (d - 1) - nullspdim
if (nlast.small.sphere > mx)
{
cat("!!! Error from pns(): \n")
cat("!!! nlast.small.sphere should be <= ",
mx,
" for this data. \n",
sep = "")
return(NULL)
}
pvalues = NaN
if (nlast.small.sphere != mx)
{
for (i in (nullspdim + 1):(d - 1 - nlast.small.sphere))
{
sp = getSubSphere(x = currentSphere, geodesic = "great")
center = sp\$center
r = sp\$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res

cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)

NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)

###########
if (nrow(currentSphere) == 3)
{
PNS\$spherePNS = t(currentSphere)
}

if (nrow(currentSphere) == 2)
{
PNS\$circlePNS = t(cur.proj)
}
#############################

}
}
if (nlast.small.sphere != 0)
{
for (i in (d - nlast.small.sphere):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center = sp\$center
r = sp\$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res

cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)

NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)

###########
if (nrow(currentSphere) == 3)
{
PNS\$spherePNS = t(currentSphere)
}

if (nrow(currentSphere) == 2)
{
PNS\$circlePNS = t(cur.proj)
}
#############################

}
}
} else {
print("!!! Error from pns():")
print("!!! sphere.type must be 'seq.test', 'small', 'great', 'BIC', or 'bi.sphere'")
print("!!!   Terminating execution     ")
return(NULL)
}
orthaxis[[d]] = meantheta
resmat[d,] = mod(S1toRadian - meantheta + pi, 2 * pi) - pi

if (output){
par(
mfrow = c(1, 1),
mar = c(4, 4, 1, 1),
mgp = c(2.5, 1, 0),
cex = 0.8
)
plot(
currentSphere[1,],
currentSphere[2,],
xlab = "",
ylab = "",
xlim = c(-1, 1),
ylim = c(-1, 1),
asp = 1
)
abline(h = 0, v = 0)
points(
cos(meantheta),
sin(meantheta),
pch = 1,
cex = 3,
col = "black",
lwd = 5
)
abline(
a = 0,
b = sin(meantheta) / cos(meantheta),
lty = 3
)
l = mod(S1toRadian - meantheta + pi, 2 * pi) - pi
points(
pch = 4,
cex = 3,
col = "blue"
)
points(
pch = 4,
cex = 3,
col = "red"
)
legend(
"topright",
legend = c("Geodesic mean", "Max (+)ve from mean", "Min (-)ve from mean"),
col = c("black", "blue", "red"),
pch = c(1, 4, 4)
)
{
cat("\n")
cat(
"length of BLUE from geodesic mean : ",
max(l),
" (",
round(max(l) * 180 / pi),
" degree)",
"\n",
sep = ""
)
cat(
"length of RED from geodesic mean : ",
min(l),
" (",
round(min(l) * 180 / pi),
" degree)",
"\n",
sep = ""
)
cat("\n")
}
}
for (i in 1:(d - 1))
{
}
resmat = flipud0(repmat(matrix(radii, ncol = 1), 1, n) * resmat)

if (d>1){
if (output){
### plot points on the 3D sphere (red), with 2D projection (blue)
rgl.sphgrid1()
sphere1.f(col="white",alpha=0.6)
}
yy <- orthaxis[[d-1]]
xx <- c(-yy[2], yy[1] , yy[3])
c1<-Enorm( c(xx[1],xx[2],xx[3])- c(-PNS\$circlePNS[1,2],PNS\$circlePNS[1,1],PNS\$circlePNS[1,3]))
costheta<- 1 - c1^2/2
angle<-(1:201)/(200)*2*pi
centre<- xx*costheta
A<- xx-centre
B<- diag(3)-A%*%t(A)/Enorm(A)**2
bv<-eigen(B)\$vectors
b1<-bv[,1]
b2<-bv[,2]
cc<- sin(acos(costheta))* ( cos(angle)%*%t(b1) + sin(angle)%*%t(b2) ) + rep(1,times=201)%*%t(centre)
if (output){
lines3d(cc,col=3,lwd=2)
}

######
if (output){
lines3d(cc,col=3,lwd=2)
sum<-0
for (i in 1:n){
sum=sum+  ( acos( cc%*%c(-PNS\$circlePNS[i,2],PNS\$circlePNS[i,1],PNS\$circlePNS[i,3])) )**2
}
mean0angle<-which.min(sum[1:200])/200*2*pi
meanpt<- sin(acos(costheta))* ( cos(mean0angle)%*%t(b1) + sin(mean0angle)%*%t(b2) ) +t(centre)
}
###########

}
PNS\$scores = t(resmat)
PNS\$pnscircle <- cbind( cbind( cc[,2],-cc[,1]) , cc[,3])
PNS\$orthaxis = orthaxis
PNS\$dist = dist
PNS\$pvalues = pvalues
PNS\$ratio = ratio
PNS\$basisu = NULL
PNS\$mean = c(PNSe2s(matrix(0, d, 1), PNS))

meanplot <- c( -PNS\$mean[2],PNS\$mean[1],PNS\$mean[3] )

if (output){
}

if (sphere.type == "seq.test")
{
PNS\$sphere.type = "seq.test"
} else if (sphere.type == "small") {
PNS\$sphere.type = "small"
} else if (sphere.type == "great") {
PNS\$sphere.type = "great"
} else if (sphere.type == "BIC") {
PNS\$sphere.type = "BIC"
} else if (sphere.type == "bi.sphere") {
PNS\$sphere.type = "bi.sphere"
}
varPNS = apply(abs(resmat) ^ 2, 1, sum) / n
total = sum(varPNS)
propPNS = varPNS / total * 100
return(list(
resmat = resmat,
PNS = PNS,
percent = propPNS
))
}

#high-res sphere plot
sphere1.f <- function(x0 = 0, y0 = 0, z0 = 0, r = 1, n = 101, ...){
f <- function(s,t){
cbind(   r * cos(t)*cos(s) + x0,
r *        sin(s) + y0,
r * sin(t)*cos(s) + z0)
}
persp3d(f, slim = c(-pi/2,pi/2), tlim = c(0, 2*pi), n = n, add = T, ...)
}

#adapted from the package "sphereplot" to remove text and axes (Aaron Robotham)
rgl.sphgrid1 <-
function (radius = 1, col.long = "red", col.lat = "blue", deggap = 15,
{
open3d()
}
for (lat in seq(-90, 90, by = deggap)) {
if (lat == 0) {
col.grid = "grey50"
}
else {
col.grid = "grey"
}
plot3d(sph2car1(long = seq(0, 360, len = 100), lat = lat,
type = "l")
}
for (long in seq(0, 360 - deggap, by = deggap)) {
if (long == 0) {
col.grid = "grey50"
}
else {
col.grid = "grey"
}
plot3d(sph2car1(long = long, lat = seq(-90, 90, len = 100),
type = "l")
}
if (longtype == "H") {
scale = 15
}
if (longtype == "D") {
scale = 1
}
# rgl.sphtext(long = 0, lat = seq(-90, 90, by = deggap), radius = radius,
#              text = seq(-90, 90, by = deggap), deg = TRUE, col = col.lat)
# rgl.sphtext(long = seq(0, 360 - deggap, by = deggap), lat = 0,
#             radius = radius, text = seq(0, 360 - deggap, by = deggap)/scale,
#              deg = TRUE, col = col.long)
#    lines3d(c(0, 0), c(0, max(radpretty)), c(0, 0), col = "grey50")
#                                                        0, radius/50), col = "grey50")
#            col = "darkgreen")
}
}
}
sph2car1<-function (long, lat, radius = 1, deg = TRUE)
{
if (is.matrix(long) || is.data.frame(long)) {
if (ncol(long) == 1) {
long = long[, 1]
}
else if (ncol(long) == 2) {
lat = long[, 2]
long = long[, 1]
}
else if (ncol(long) == 3) {
lat = long[, 2]
long = long[, 1]
}
}
if (missing(long) | missing(lat)) {
stop("Missing full spherical 3D input data.")
}
if (deg) {
long = long * pi/180
lat = lat * pi/180
}
return = cbind(x = radius * cos(long) * cos(lat), y = radius *
sin(long) * cos(lat), z = radius * sin(lat))
}

#==================================================================================
pns4pc = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 1,
n.pc = 2)
{
if (n.pc < 2)
{
stop("Error: n.pc should be >= 2.")
}
out = pc2sphere2(x = x, n.pc = n.pc)
spheredata = t(out\$spheredata)
GPAout = out\$GPAout
pns.out = pns(
x = spheredata,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere
)
pns.out\$percent = pns.out\$percent * sum(GPAout\$percent[1:n.pc]) / 100
pns.out\$GPAout = GPAout
pns.out\$spheredata = spheredata
return(pns.out)
}

pns.pc = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 0,
n.pc = 0)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n.pc == 0)
{
GPAout = procGPA(
x = x,
scale = TRUE,
reflect = FALSE,
tangentcoords = "partial",
distances = FALSE
)
spheredata = matrix(NA, k * m, n)
for (i in 1:n)
{
spheredata[, i] = c(GPAout\$rotated[, , i])
}
pns.out = pns(
x = spheredata,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere
)
resmat = pns.out\$resmat
PNS = pns.out\$PNS
npts = 200
prinarc1 = get.prinarc(
resmat,
PNS,
arc = 1,
n = npts,
boundary.data = FALSE
)
prinarc2 = get.prinarc(
resmat,
PNS,
arc = 2,
n = npts,
boundary.data = FALSE
)
prinarc1.ar = array(NA, c(k, m, npts))
prinarc2.ar = array(NA, c(k, m, npts))
for (i in 1:npts)
{
prinarc1.ar[, , i] = matrix(prinarc1[, i], nrow = k)
prinarc2.ar[, , i] = matrix(prinarc2[, i], nrow = k)
}
scores.prinarc1 = shape.pcscores.partial(PCAout = GPAout, x = prinarc1.ar)
scores.prinarc2 = shape.pcscores.partial(PCAout = GPAout, x = prinarc2.ar)
out = pns.out
out\$GPAout = GPAout
out\$scores.prinarc1 = scores.prinarc1
out\$scores.prinarc2 = scores.prinarc2
} else {
pns.out = pns4pc(
x = x,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere,
n.pc = n.pc
)
GPAout = pns.out\$GPAout
resmat = pns.out\$resmat
PNS = pns.out\$PNS
npts = 200
prinarc1 = get.prinarc(
resmat,
PNS,
arc = 1,
n = npts,
boundary.data = FALSE
)
prinarc2 = get.prinarc(
resmat,
PNS,
arc = 2,
n = npts,
boundary.data = FALSE
)
scores.prinarc1 = matrix(NA, npts, n.pc)
scores.prinarc2 = matrix(NA, npts, n.pc)
for (g in 1:npts)
{
size1 = acos(prinarc1[1, g])
size2 = acos(prinarc2[1, g])
scores.prinarc1[g,] = prinarc1[2:(n.pc + 1), g] / (sin(size1) / size1)
scores.prinarc2[g,] = prinarc2[2:(n.pc + 1), g] / (sin(size2) / size2)
}
out = pns.out
out\$scores.prinarc1 = scores.prinarc1
out\$scores.prinarc2 = scores.prinarc2
}
return(out)
}

#==================================================================================
rotMat = function(b, a = NULL, alpha = NULL)
{
if (is.matrix(b))
{
if (min(dim(b)) == 1)
{
b = c(b)
} else {
stop("Error: b should be a unit vector.")
}
}
d = length(b)
b = b / norm(b, type = "2")

if (is.null(a) & is.null(alpha))
{
a = c(rep(0, d - 1), 1)
alpha = acos(sum(a * b))
} else if (!is.null(a) & is.null(alpha)) {
alpha = acos(sum(a * b))
} else if (is.null(a) & !is.null(alpha)) {
a = c(rep(0, d - 1), 1)
}
if (abs(sum(a * b) - 1) < 1e-15)
{
rot = diag(d)
return(rot)
}
if (abs(sum(a * b) + 1) < 1e-15)
{
rot = -diag(d)
return(rot)
}
c = b - a * sum(a * b)
c = c / norm(c, type = "2")
A = a %*% t(c) - c %*% t(a)
rot = diag(d) + sin(alpha) * A + (cos(alpha) - 1) * (a %*% t(a) + c %*% t(c))
return(rot)
}

#==================================================================================
ExpNPd = function(x)
{
if (is.vector(x))
{
x = as.matrix(x)
}
d = nrow(x)
nv = sqrt(apply(x ^ 2, 2, sum))
Exppx = rbind(matrix(rep(sin(nv) / nv, d), nrow = d, byrow = T) * x, cos(nv))
Exppx[, nv < 1e-16] = repmat(matrix(c(rep(0, d), 1)), 1, sum(nv < 1e-16))
return(Exppx)
}

#==================================================================================
LogNPd = function(x)
{
n = ncol(x)
d = nrow(x)
scale = acos(x[d,]) / sqrt(1 - x[d,] ^ 2)
scale[is.nan(scale)] = 1
Logpx = repmat(t(scale), d - 1, 1) * x[-d,]
return(Logpx)
}

#==================================================================================
objfn = function(center, r, x)
{
return(mean((acos(t(
center
) %*% x) - r) ^ 2))
}

#==================================================================================
getSubSphere = function(x, geodesic = "small")
{
svd.x = svd(x)
initialCenter = svd.x\$u[, ncol(svd.x\$u)]
c0 = initialCenter
TOL = 1e-10
cnt = 0
err = 1
n = ncol(x)
d = nrow(x)
Gnow = 1e+10
while (err > TOL)
{
c0 = c0 / norm(c0, type = "2")
rot = rotMat(c0)
TpX = LogNPd(rot %*% x)
fit = sphereFit(
x = TpX,
initialCenter = rep(0, d - 1),
geodesic = geodesic
)
newCenterTp = fit\$center
r = fit\$r
if (r > pi)
{
r = pi / 2
svd.TpX = svd(TpX)
newCenterTp = svd.TpX\$u[, ncol(svd.TpX\$u)] * pi / 2
}
newCenter = ExpNPd(newCenterTp)
center = solve(rot, newCenter)
Gnext = objfn(center, r, x)
err = abs(Gnow - Gnext)
Gnow = Gnext
c0 = center
cnt = cnt + 1
if (cnt > 30)
{
break
}
}
i1save = list()
i1save\$Gnow = Gnow
i1save\$center = center
i1save\$r = r
initialCenter = U[, ncol(U)]
c0 = initialCenter
TOL = 1e-10
cnt = 0
err = 1
n = ncol(x)
d = nrow(x)
Gnow = 1e+10

while (err > TOL)
{
c0 = c0 / norm(c0, type = "2")
rot = rotMat(c0)
TpX = LogNPd(rot %*% x)
fit = sphereFit(
x = TpX,
initialCenter = rep(0, d - 1),
geodesic = geodesic
)
newCenterTp = fit\$center
r = fit\$r
if (r > pi)
{
r = pi / 2
svd.TpX = svd(TpX)
newCenterTp = svd.TpX\$u[, ncol(svd.TpX\$u)] * pi / 2
}
newCenter = ExpNPd(newCenterTp)
center = solve(rot, newCenter)
Gnext = objfn(center, r, x)
err = abs(Gnow - Gnext)
Gnow = Gnext
c0 = center
cnt = cnt + 1
if (cnt > 30)
{
break
}
}
if (i1save\$Gnow == min(Gnow, i1save\$Gnow))
{
center = i1save\$center
r = i1save\$r
}
if (r > pi / 2)
{
center = -center
r = pi - r
}
return(list(center = c(center), r = r))
}

#==================================================================================
LRTpval = function(resGREAT, resSMALL, n)
{
chi2 = max(n * log(sum(resGREAT ^ 2) / sum(resSMALL ^ 2)), 0)
return(pchisq(
q = chi2,
df = 1,
lower.tail = FALSE
))
}

#==================================================================================
vMFtest = function(x, R = 100)
{
d = nrow(x)
n = ncol(x)
sumx = apply(x, 1, sum)
rbar = norm(sumx, "2") / n
muMLE = sumx / norm(sumx, "2")
kappaMLE = (rbar * d - rbar ^ 3) / (1 - rbar ^ 2)
sp = getSubSphere(x = x, geodesic = "small")
center.s = sp\$center
r.s = sp\$r
xi_vec = rep(0, R)
for (r in 1:R)
{
rdata = randvonMisesFisherm(d, n, kappaMLE)
sp = getSubSphere(x = rdata, geodesic = "small")
center.s = sp\$center
r.s = sp\$r
}
pvalue = mean(xi_vec > xi_sample)
return(pvalue)
}

#==================================================================================
geodmeanS1 = function(theta)
{
n = length(theta)
meancandi = mod(mean(theta) + 2 * pi * (0:(n - 1)) / n, 2 * pi)
theta = mod(theta, 2 * pi)
geodvar = rep(0, n)
for (i in 1:n)
{
v = meancandi[i]
dist2 = apply(cbind((theta - v) ^ 2, (theta - v + 2 * pi) ^ 2, (v - theta + 2 * pi) ^
2), 1, min)
geodvar[i] = sum(dist2)
}
m = min(geodvar)
ind = which.min(geodvar)
geodmean = mod(meancandi[ind], 2 * pi)
geodvar = geodvar[ind] / n
return(list(geodmean = geodmean, geodvar = geodvar))
}

#==================================================================================
PNSe2s = function(resmat, PNS)
{
dm = nrow(resmat)
n = ncol(resmat)
NSOrthaxis = rev(PNS\$orthaxis[1:(dm - 1)])
NSradius = flipud0(matrix(PNS\$dist, ncol = 1))
geodmean = PNS\$orthaxis[[dm]]
res = resmat / repmat(flipud0(matrix(PNS\$radii, ncol = 1)), 1, n)
T = t(rotMat(NSOrthaxis[[1]])) %*%
rbind(repmat(sin(NSradius[1] + matrix(res[2,], nrow = 1)), 2, 1) *
rbind(cos(geodmean + res[1,]), sin(geodmean + res[1,])),
if (dm > 2)
{
for (i in 1:(dm - 2))
{
T = t(rotMat(NSOrthaxis[[i + 1]])) %*%
res[i + 2,], nrow = 1
)), 2 + i, 1) * T,
cos(NSradius[i + 1] + res[i + 2,]))
}
}
if (!is.null(PNS\$basisu))
{
T = PNS\$basisu %*% T
}
return(T)
}

#==================================================================================
PNSs2e = function(spheredata, PNS)
{
if (nrow(spheredata) != length(PNS\$mean))
{
cat(" Error from PNSs2e() \n")
cat(" Dimensions of the sphere and PNS decomposition do not match")
return(NULL)
}
if (!is.null(PNS\$basisu))
{
spheredata = t(PNS\$basisu) %*% spheredata
}
kk = nrow(spheredata)
n = ncol(spheredata)
Res = matrix(0, kk - 1, n)
currentSphere = spheredata
for (i in 1:(kk - 2))
{
v = PNS\$orthaxis[[i]]
r = PNS\$dist[i]
res = acos(t(v) %*% currentSphere) - r
Res[i,] = res
NestedSphere = rotMat(v) %*% currentSphere
currentSphere = as.matrix(NestedSphere[1:(kk - i),]) /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), kk - i, 1)
}
devS1 = mod(S1toRadian - rev(PNS\$orthaxis)[[1]] + pi, 2 * pi) - pi
Res[kk - 1,] = devS1
EuclidData = flipud0(repmat(PNS\$radii, 1, n) * Res)
return(EuclidData)
}

#==================================================================================
randvonMisesFisherm = function(m, n, kappa, mu = NULL)
{
if (is.null(mu))
{
muflag = FALSE
} else {
muflag = TRUE
}
if (m < 2)
{
print("Message from randvonMisesFisherm(): dimension m must be > 2")
print("Message from randvonMisesFisherm(): Set m to be 2")
m = 2
}
if (kappa < 0)
{
print("Message from randvonMisesFisherm(): kappa must be >= 0")
print("Message from randvonMisesFisherm(): Set kappa to be 0")
kappa = 0
}
b = (-2 * kappa + sqrt(4 * kappa ^ 2 + (m - 1) ^ 2)) / (m - 1)
x0 = (1 - b) / (1 + b)
c = kappa * x0 + (m - 1) * log(1 - x0 ^ 2)
nnow = n
w = c()
while (TRUE)
{
ntrial = max(round(nnow * 1.2), nnow + 10)
Z = rbeta(n = ntrial,
shape1 = (m - 1) / 2,
shape2 = (m - 1) / 2)
U = runif(ntrial)
W = (1 - (1 + b) * Z) / (1 - (1 - b) * Z)

indicator = kappa * W + (m - 1) * log(1 - x0 * W) - c >= log(U)
if (sum(indicator) >= nnow)
{
w1 = W[indicator]
w = c(w, w1[1:nnow])
break
} else {
w = c(w, W[indicator])
nnow = nnow - sum(indicator)
}
}
V = UNIFORMdirections(m - 1, n)
X = rbind(repmat(sqrt(1 - matrix(w, nrow = 1) ^ 2), m - 1, 1) * V, matrix(w, nrow = 1))
if (muflag)
{
mu = mu / norm(mu, "2")
X = t(rotMat(mu)) %*% X
}
return(X)
}

#==================================================================================
UNIFORMdirections = function(m, n)
{
V = matrix(0, m, n)
nr = matrix(rnorm(m * n), nrow = m)
for (i in 1:n)
{
while (TRUE)
{
ni = sum(nr[, i] ^ 2)
if (ni < 1e-10)
{
nr[, i] = rnorm(m)
} else {
V[, i] = nr[, i] / sqrt(ni)
break
}
}
}
return(V)
}

#==================================================================================
trans.subsphere = function(x, center)
{
return(repmat(1 / sqrt(1 - (t(
center
) %*% x) ^ 2), length(center) - 1, 1) *
(rotMat(center)[-length(center),] %*% x))
}

#==================================================================================
get.prinarc.value = function(PNS, arc, res)
{
d = length(PNS\$orthaxis)
n = length(res)
prinarc = matrix(NA, d + 1, n)
for (g in 1:n)
{
newres = matrix(0, d, 1)
newres[arc] = res[g]
T = PNSe2s(newres, PNS)
prinarc[, g] = T
}
return(prinarc)
}

#==================================================================================
get.prinarc = function(resmat, PNS, arc, n, boundary.data = FALSE)
{
d = nrow(resmat)
if (boundary.data)
{
mn = min(resmat[arc,])
mx = max(resmat[arc,])
} else {
mn = -pi * tail(PNS\$radii, arc)[1]
mx = pi * tail(PNS\$radii, arc)[1]
}
prinarcgrid = seq(mn, mx, length = n)
prinarc = matrix(NA, d + 1, n)
for (g in 1:n)
{
newres = matrix(0, d, 1)
newres[arc] = prinarcgrid[g]
T = PNSe2s(newres, PNS)
prinarc[, g] = T
}
return(prinarc)
}

#==================================================================================
get.prinarc.subsphere = function(resmat,
PNS,
arc,
n,
subsphere = arc,
boundary.data = FALSE)
{
if (subsphere < arc)
{
stop("Error: subsphere >= arc.")
}
if (subsphere < 1)
{
stop("Error: subsphere >= 1.")
}
prinarc = get.prinarc(
resmat = resmat,
PNS = PNS,
arc = arc,
n = n,
boundary.data = boundary.data
)
d = nrow(resmat)
prinarc.sub = prinarc
if (subsphere < d)
{
for (i in 1:(d - subsphere))
{
prinarc.sub = trans.subsphere(x = prinarc.sub, center = PNS\$orthaxis[[i]])
}
}
return(prinarc.sub)
}

#==================================================================================
get.data.subsphere = function(resmat, PNS, x, subsphere)
{
if (subsphere < 1)
{
stop("Error: subsphere >= 1.")
}
d = nrow(resmat)
x.sub = x
if (subsphere < d)
{
for (i in 1:(d - subsphere))
{
x.sub = trans.subsphere(x = x.sub, center = PNS\$orthaxis[[i]])
}
}
return(x.sub)
}

#==================================================================================
mod = function(x, y)
{
return(x %% y)
}

#==================================================================================
repmat = function(x, m, n)
{
return(kronecker(matrix(1, m, n), x))
}

#==================================================================================
flipud0 = function(x)
{
return(apply(x, 2, rev))
}

#==================================================================================
sphere.obj = function(center, x, is.greatCircle)
{
di = sqrt(apply((x - repmat(
matrix(center, ncol = 1), 1, ncol(x)
)) ^ 2, 2, sum))
if (is.greatCircle)
{
r = pi / 2
} else {
r = mean(di)
}
sum((di - r) ^ 2)
}

#==================================================================================
sphere.res = function(center, x, is.greatCircle)
{
center = c(center)
xmc = x - center
di = sqrt(apply(xmc ^ 2, 2, sum))
if (is.greatCircle)
{
r = pi / 2
} else {
r = mean(di)
}
(di - r)
}

#==================================================================================
sphere.jac = function(center, x, is.greatCircle)
{
center = c(center)
xmc = x - center
di = sqrt(apply(xmc ^ 2, 2, sum))
di.vj = -xmc / repmat(matrix(di, nrow = 1), length(center), 1)
if (is.greatCircle)
{
c(t(di.vj))
} else {
r.vj = apply(di.vj, 1, mean)
c(t(di.vj - repmat(matrix(r.vj, ncol = 1), 1, ncol(x))))
}
}

#==================================================================================
sphereFit = function(x,
initialCenter = NULL,
geodesic = "small")
{
if (is.null(initialCenter))
{
initialCenter = apply(x, 1, mean)
}
op = nls.lm(
par = initialCenter,
fn = sphere.res,
jac = sphere.jac,
x = x,
is.greatCircle = ifelse(geodesic == "great", TRUE, FALSE),
control = nls.lm.control(maxiter = 1000)
)
center = coef(op)
di = sqrt(apply((x - repmat(
matrix(center, ncol = 1), 1, ncol(x)
)) ^ 2, 2, sum))
if (geodesic == "great")
{
r = pi / 2
} else {
r = mean(di)
}

list(center = center, r = r)
}

#==================================================================================
tr = function(x)
{
return(sum(diag(x)))
}

#==================================================================================
Enormalize = function(x)
{
return(x / Enorm(x))
}

#==================================================================================
sphere2pcscore = function(x)
{
n = nrow(x)
p = ncol(x)
scores = matrix(NA, n, p - 1)
for (i in 1:n)
{
size = acos(x[i, 1])
scores[i,] = (size / sin(size)) * x[i, 2:p]
}
return(scores)
}

#==================================================================================
pcscore2sphere = function(n.pc, X.hat, S, Tan, V)
{
d = nrow(Tan)
n = ncol(Tan)
W = matrix(NA, d, n)
for (i in 1:n)
{
W[, i] = acos(tr(S[, , i] %*% t(X.hat))) * Tan[, i] / sqrt(sum(Tan[, i] ^
2))
}
lambda = matrix(NA, n, d)
for (i in 1:n)
{
for (j in 1:d)
{
lambda[i, j] = sum(W[, i] * V[, j])
}
}
U = matrix(0, n, d)
for (i in 1:n)
{
for (j in 1:n.pc)
{
U[i,] = U[i,] + lambda[i, j] * V[, j]
}
}
S.star = matrix(NA, n, n.pc + 1)
for (i in 1:n)
{
U.norm = sqrt(sum(U[i,] ^ 2))
S.star[i,] = c(cos(U.norm),
sin(U.norm) / U.norm * lambda[i, 1:n.pc])
}
return(S.star)
}

pcscore2sphere2 = function(n.pc, X.hat, S, Tan, V)
{
d = nrow(Tan)
n = ncol(Tan)
W = matrix(NA, d, n)

if (n.pc > min(d,n) )
{
stop("Error: n.pc must be <= min(n,d)")
}

for (i in 1:n)
{
W[, i] = acos(tr(S[, , i] %*% t(X.hat))) * Tan[, i] / sqrt(sum(Tan[, i] ^
2))
}
lambda = matrix(NA, n, d)
for (i in 1:n)
{
for (j in 1:n.pc)
{
lambda[i, j] = sum(W[, i] * V[, j])
}
}
U = matrix(0, n, d)
for (i in 1:n)
{
for (j in 1:n.pc)
{
U[i,] = U[i,] + lambda[i, j] * V[, j]
}
}
S.star = matrix(NA, n, n.pc + 1)
for (i in 1:n)
{
U.norm = sqrt(sum(U[i,] ^ 2))
S.star[i,] = c(cos(U.norm),
sin(U.norm) / U.norm * lambda[i, 1:n.pc])
}
return(S.star)
}

#==================================================================================
pc2sphere = function(x, n.pc)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n.pc < ((k - 1) * m))
{
stop("Error: n.pc must be >= (k - 1) * m.")
}
GPAout = procGPA(
x = x,
scale = TRUE,
reflect = FALSE,
tangentcoords = "partial",
distances = FALSE
)
cat(
"First ",
n.pc,
" principal components explain ",
round(sum(GPAout\$percent[1:n.pc])),
"% of total variance. \n",
sep = ""
)
H = defh(k - 1)
X.hat = H %*% GPAout\$mshape
S = array(NA, c(k - 1, m, n))
for (i in 1:n)
{
S[, , i] = H %*% GPAout\$rotated[, , i]
}
T.c = GPAout\$tan - apply(GPAout\$tan, 1, mean)
out = pcscore2sphere(
n.pc = n.pc,
X.hat = X.hat,
S = S,
Tan = T.c,
V = GPAout\$pcar
)
return(list(spheredata = out, GPAout = GPAout))
}

#==================================================================================
rot.mat = function(Y,
X,
reflect = FALSE,
center = TRUE)
{
svd.out = svd(t(X) %*% Y)
R = svd.out\$u %*% t(svd.out\$v)
if (!reflect)
{
if (det(R) < 0)
{
u = svd.out\$u
v = svd.out\$v
if (det(u) < 0)
{
u[, dim(u)[2]] = -u[, dim(u)[2]]
} else if (det(v) < 0) {
v[, dim(v)[2]] = -v[, dim(v)[2]]
}
R = u %*% t(v)
}
}
return(R)
}

#==================================================================================
Procrustes.dist.full = function(x1, x2)
{
m = ncol(x1)
z1 = preshape(x1)
z2 = preshape(x2)
Q = t(z1) %*% z2 %*% t(z2) %*% z1
ev = eigen(Q)\$values
sign = ifelse(det(t(z1) %*% z2) >= 0, 1,-1)
dF = sqrt(abs(1 - sum(sqrt(abs(
ev[1:(m - 1)]
)), sign * sqrt(abs(
ev[m]
))) ^ 2))
R = rot.mat(
Y = z2,
X = z1,
reflect = FALSE,
center = FALSE
)
scale = sum(svd(t(z1) %*% z2)\$d)
return(list(dF = dF, R = R, scale = scale))
}

#==================================================================================
tangent.coords.partial = function(x, p)
{
k = nrow(x)
m = ncol(x)
if (abs(norm(p, "F") - 1) > 1e-15)
{
print("||p|| is not 1. Normalised one is used.")
p = Enormalize(p)
}
tmp = Procrustes.dist.full(x, p)
R = tmp\$R
scale = tmp\$scale
pre.p = preshape(p)
pre.x = preshape(x)
ident = diag(k * m - m)
tan = (ident - matrix(pre.p) %*% t(c(pre.p))) %*% c(pre.x %*% R)
tan.scale = (ident - matrix(pre.p) %*% t(c(pre.p))) %*% c(pre.x %*% R * scale)
return(list(
tan = c(tan),
tan.scale = c(tan.scale),
R = R,
scale = scale
))
}

#==================================================================================
shape.pcscores = function(PCAout, x, tangentcoords = "partial")
{
if (tangentcoords == "partial")
{
if (abs(norm(PCAout\$mshape, "F") - 1) > 1e-15)
{
print("||PCAout\$mshape|| is not 1. Normalised one is used.")
mshape = Enormalize(PCAout\$mshape)
} else {
mshape = PCAout\$mshape
}
if (abs(norm(x, "F") - 1) > 1e-15)
{
print("||x|| is not 1. Normalised one is used.")
x = Enormalize(x)
}
opa.out = procOPA(mshape, x, scale = FALSE)
matched = opa.out\$Bhat
tan.out = tangent.coords.partial(matched, mshape)
mean.tan = apply(PCAout\$tan, 1, mean)
scores = t(tan.out\$tan - mean.tan) %*% PCAout\$pcar
scores.scale = t(tan.out\$tan.scale - mean.tan) %*% PCAout\$pcar
return(
list(
rotated = matched,
tan = tan.out\$tan,
tan.scale = tan.out\$tan.scale,
scores = c(scores),
scores.scale = c(scores.scale)
)
)
}
}

#==================================================================================
shape.pcscores.partial = function(PCAout, x)
{
n = dim(x)[3]
scores = c()
for (i in 1:n)
{
s = shape.pcscores(PCAout, x[, , i], tangentcoords = "partial")
scores = rbind(scores, s\$scores)
}
return(scores)
}

#==================================================================================
plotshapes3d.pns = function(x,
type = "p",
col = "black",
size = 5,
aspect = "iso",
joinline = TRUE,
col.joinline = "#d4d2d2",
lwd.joinline = 0.5,
tick = FALSE,
labels.tick = FALSE,
xlab = "",
ylab = "",
zlab = "")
{
k = dim(x)[1]
n = dim(x)[3]
aa = c()
bb = c()
cc = c()
for (i in 1:n)
{
aa = c(aa, x[, 1, i])
bb = c(bb, x[, 2, i])
cc = c(cc, x[, 3, i])
}
xlim = range(aa)
ylim = range(bb)
zlim = range(cc)
plot3d(
x[, , 1],
type = "n",
xlab = "",
ylab = "",
zlab = "",
box = FALSE,
axes = FALSE,
aspect = aspect,
xlim = xlim,
ylim = ylim,
zlim = zlim
)
for (i in 1:n)
{
plot3d(
x[, , i],
type = type,
col = col,
size = size,
)
}
if (tick)
{
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
} else {

}
r = cbind(xlim, ylim, zlim)
pos = r[2,] + apply(r, 2, diff) / 20
text3d(pos[1], 0, 0, texts = xlab, cex = 0.8)
text3d(0, pos[2], 0, texts = ylab, cex = 0.8)
text3d(0, 0, pos[3], texts = zlab, cex = 0.8)
if (joinline)
{
for (i in 1:n)
{
lines3d(x[, , i], col = col.joinline, lwd = lwd.joinline)
}
}
}

#==================================================================================
Plot3D = function(x,
type = "s",
col = "black",
size = 1.2,
aspect = "iso",
joinline = FALSE,
col.joinline = "#d4d2d2",
lwd.joinline = 0.5,
tick = TRUE,
tick.boundary = FALSE,
labels.tick = TRUE,
xlab = "",
ylab = "",
zlab = "")
{
n = nrow(x)
plot3d(
x,
type = "n",
xlab = "",
ylab = "",
zlab = "",
box = FALSE,
axes = FALSE,
aspect = aspect
)
plot3d(
x,
type = type,
col = col,
size = size,
)
if (tick)
{
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
}
if (tick.boundary)
{
tks = pretty(x[, 1], n = 10)
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
tks = pretty(x[, 2], n = 10)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
tks = pretty(x[, 3], n = 10)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
}
r = apply(x, 2, range)
pos = r[2,] + apply(r, 2, diff) / 20
text3d(pos[1], 0, 0, texts = xlab, cex = 0.8)
text3d(0, pos[2], 0, texts = ylab, cex = 0.8)
text3d(0, 0, pos[3], texts = zlab, cex = 0.8)
if (joinline)
{
lines3d(x, col = col.joinline, lwd = lwd.joinline)
}
}

#==================================================================================
col2RGB = function(col, alpha = 255)
{
n = length(col)
out = c()
for (i in 1:n)
{
out[i] = rgb(
red = col2rgb(col[i])[1],
green = col2rgb(col[i])[2],
blue = col2rgb(col[i])[3],
alpha = alpha,
maxColorValue = 255
)
}
return(out)
}

#==================================================================================
project.subsphere = function(x, center, r)
{
n = ncol(x)
d = nrow(x)
x.proj = matrix(NA, d, n)
for (i in 1:n)
{
rho = acos(sum(x[, i] * center))
x.proj[, i] = (sin(r) * x[, i] + sin(rho - r) * center) / sin(rho)
}
return(x.proj)
}

##############################################################end of PNS###########

##### Penalised Euclidean Distance Regression

#==================================================================================
ped <- function(X, Y, method = c("AIC")) {
if (method == "AIC") {
aicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.35, 1.5)) {
out   <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (out\$aic < aicmin) {
minout <- out
mincofp <- cofp
aicmin <- out\$aic
}
}
}
out <- minout
}

if (method == "BIC") {
bicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.35, 1.5)) {
out   <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (out\$bic < bicmin) {
minout <- out
mincofp <- cofp
bicmin <- out\$bic
}
}
}
out <- minout
}

if (method == "khat") {
aicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.25, 1.5)) {
out   <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (-out\$khat < aicmin) {
minout <- out
mincofp <- cofp
aicmin <- -out\$khat
}
}
}
out <- minout
}
if (method == "CV") {
n <- length(Y)
cvmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.25, 1.5)) {
cverr <- 0
for (jj in 1:10) {
subsample <- ((jj - 1) * 10 + 1):((jj - 1) * 10 + 10)
out   <-
pedreg(
X[-subsample, ],
Y[-subsample],
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
cverr <-
cverr + Enorm(Y[subsample] - out\$intercept - X[subsample, ] %*% out\$betahat) **
2
}
if (cverr < cvmin) {
minout <- out
minlam <- lam
mincofp <- cofp
cvmin <- cverr
}
}
}
out <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = mincofp,
lambdainit = minlam
)
}

out1 <- list(
betahat = 0,
yhat = 0,
lambda = 0,
coef = 0,
resid = 0
)
out1\$intercept <- out\$intercept
out1\$coef <- c(out\$intercept, out\$betahat)
out1\$betahat <- out\$betahat
out1\$lambda <- out\$lambda
out1\$delta <- mincofp
out1\$yhat <- out\$yhat
out1\$resid <- Y - out\$yhat

out1
}

###########################function for PED#####################

#==================================================================================
pedreg <-
function(X,
Y,
constc0 = 1.1,
constc1 = 1.35,
alpha = 0.05,
LMM = 50,
MIT = 10000,
NUM_METHOD = 1,
nlambda = 1,
lambdamax = 1,
PLOT = TRUE,
BIC = FALSE,
lambdainit = 1) {
# NUM_METHOD = 1 = L-BFGS-B
# LMM = Parameter M in L-BFGS method 1
# MIT = Max iterations for optimization

p <- dim(X)[2]
n <- dim(X)[1]
constc <- constc0
Ymean <- mean(Y)
Ysd <- sd(Y)
Yinit <- Y

pinit <- p
ans0 <- rep(0, times = pinit)

Xorig <- X
Yorig <- Y

vm <- rep(0, times = ncol(X))
vsd <- rep(0, times = ncol(X))

for (i in 1:ncol(X)) {
vm[i] <- mean(X[, i])
}
for (i in 1:ncol(X)) {
vsd[i] <- sd(X[, i])
}

#standardize to sphere

X <- scale(X) / sqrt(n - 1)
Y <- scale(Y) / sqrt(n - 1)

X0 <- X
Y0 <- Y

lambdainit1 <- constc / sqrt(n - 1) * sqrt(sqrt(p)) / n * qnorm(1 - alpha /
(2 * p))
if (nlambda == 1) {
lambdainit1 <- lambdainit
}

METHOD1 <- "L-BFGS-B"
xi <- -9999999
c1 <- 1

nlam <- nlambda
betamat <- matrix(0, p, nlam)
betamat.sparse <- betamat
lambdamat <- rep(0, times = nlam)
ximat <- rep(0, times = nlam)
aic <- rep(0, times = nlam)
bic <- rep(0, times = nlam)
npar <- rep(0, times = nlam)
selectmat <- betamat
#cat(c("Lambda iteration (out of ",nlam,"):"))

for (ilam in (nlam:1)) {
#cat(c(ilam," "))

if (nlam == 1) {
lambda <- lambdainit1
}
if (nlam > 1) {
c1 <- sqrt(n) + (ilam - 1) / (nlam - 1) * 1 / lambdainit1
c1 <-
sqrt(n) + ((ilam - 1) / (nlam - 1)) * 1 / lambdainit1 * (lambdamax - lambdainit1 *
sqrt(n))
lambda <- lambdainit1 * c1
}

if (ilam == nlam) {
x0 <- rep(1 / sqrt(p), times = p)
}

if (ilam != nlam) {
x0 <- betahat + rnorm(p) / sqrt(p)
}
#x0<-rnorm(p)/sqrt(p)

pedfun <- function(pars, Y = 0, X = 0) {
p <- length(pars)
pars <- matrix(pars, p, 1)
ped <- Enorm(Y - X %*% pars) + lambda * sqrt(Enorm(pars) * sum(abs(pars)))
ped
}

pedgrad <- function(pars, X = 0, Y = 0) {
GM <- sqrt(Enorm(pars) * sum(abs(pars)))
gradL <- rep(0, times = p)

gradL <-  -t(X) %*% (Y - X %*% pars) / Enorm(Y - X %*% pars)
lambda / 2 * pars / Enorm(pars) * sum(abs(pars)) / GM +
lambda / 2 * sign(pars) * Enorm(pars) / GM ,
p,
1
)

}

if (NUM_METHOD == 1) {
repeat {
#,ndeps=1e-3,factr=1e-5,pgtol=1e-5
res2 <-
optim(
par = x0,
fn = pedfun,
method = METHOD1,
control = list(lmm = LMM, maxit = MIT),
X = X,
Y = Y
)
betahat <- res2\$par
if (res2\$convergence == 0) {
break
}
x0 <- rnorm(p) / sqrt(p)
}
}

oldxi <- xi
xi <-
sqrt(Enorm(betahat) / sum(abs(betahat))) - sqrt(n) / (constc * c1 * p ^
(1 / 4))

dif <- (xi - oldxi)

REGC <- 0.0001

betamat[, ilam] <- betahat / (Enorm(betahat) + REGC)

lambdamat[ilam] <- lambda
ximat[ilam] <- xi

ximat[ilam] <- sqrt(Enorm(betahat) / sum(abs(betahat)))

MM <- constc1 / (sqrt(n))

select <- (abs(betahat) / (Enorm(betahat) + REGC) > MM)
selectmat[, ilam] <- select
betamat.sparse[, ilam] <- betamat[, ilam]
betamat.sparse[select == FALSE, ilam] <- 0 * betamat[select == FALSE, ilam]

pp <- sum(select)
npar[ilam] <- pp
bic[ilam] <- log(Enorm(Y) ** 2 / n) * (n) + log(n) * (1)
aic[ilam] <- log(Enorm(Y) ** 2 / n) * (n) + 2 * (1)

if (sum(select) > 0) {
aa <- lm(Y ~ X[, c(1:p)[select]] - 1)
pred <- predict(aa)

# Use AIC with finite sample correction
aic[ilam] <-
log(Enorm(Y - pred) ** 2 / n) * (n)  + 2 * (pp + 1) + 2 * (pp + 1) * (pp +
2) / (n - pp - 2)
# Use AIC/BIC
bic[ilam] <- log(Enorm(Y - pred) ** 2 / n) * (n) + log(n) * (pp + 1)
aic[ilam] <- log(Enorm(Y - pred) ** 2 / n) * (n) + 2 * (pp + 1)
}

}

best <- 1
if (nlam > 1) {
###########################################choose best via AIC

best <- c(1:nlam)[aic == min(aic)][1]
select <- as.logical(selectmat[, best])
lambdaaic <- lambdamat[best]

###########################################choose best via Corollary 1 with khat
best2 <- nlam
if ((sum(ximat > 0.25)) > 0) {
best2 <- c(1:nlam)[(ximat) > 0.25][1]
}

#################### biggest xi from Corollary 1

xism <- (ximat)
if (sum(diff(xism) < 0.01) > 0) {
best2 <- c(2:nlam)[diff(xism) < 0.01][1] - 1
}

############## biggest sqrt( Enorm(beta) / norm(beta)_1 )

best2 <- c(1:nlam)[ximat == max(ximat)]

selectcor <- as.logical(selectmat[, best2])
lambdacor1 <- lambdamat[best2]

if (BIC == FALSE) {
best <- best2
select <- selectcor
}
}

################last part######estimate with reduced p###########

if (sum(select) > 0) {
X <- as.matrix(X[, select])
p <- sum(select)
p14 <- sqrt(sqrt(p))
final <- c(1:pinit)[select]
}

#######

lambda <- constc / sqrt(sqrt(p)) / sqrt(n) * qnorm(1 - alpha / (2 * p))

if (sum(select) > 0) {
x0 <- rep(1 / p, times = p)

if (NUM_METHOD == 1) {
repeat {
res2 <-
optim(
par = x0,
fn = pedfun,
method = METHOD1,
control = list(lmm = LMM, maxit = MIT),
X = X,
Y = Y
)
betahat <- res2\$par

if (res2\$convergence == 0) {
break
}

x0 <- rnorm(p) / p
}
}

ans0[final] <- betahat

}

#ind<-which(abs(ans0/Enorm(ans0))<10^(-5))
#ans0[ind]<-0

out <- list(betahat = 0,
yhat = 0,
lambda = 0)

out\$betahatscale <- ans0
out\$yhatscale <- c(X0 %*% ans0)

if (nlam > 1) {
out\$lambdacor1 <- lambdacor1
out\$lambdaaic <- lambdaaic
out\$betamat.sparse <- betamat.sparse
out\$betamat.rescale <- betamat
out\$betamat <- betamat
for (i in 1:nlam) {
out\$betamat.rescale[, i] <- c(out\$betamat[, i] / vsd) * sd(Yorig)
}
out\$lambdamat <- lambdamat
out\$ximat <- ximat
out\$MM <- MM
out\$fmax <- res2\$value
out\$npar <- npar
out\$selectmat <- selectmat
}

#use AIC
out\$aic <- aic
#use BIC
out\$bic <- bic

#use khat
out\$khat <- ximat

out\$lambdath3 <- lambdainit1 * sqrt(n)

out\$lambda <- out\$lambdacor1
if (BIC == TRUE) {
out\$lambda <- out\$lambdaaic
}

if (nlam == 1) {
out\$lambda <- lambdainit1
out\$constc1 <- constc1
}

sol <- sd(Yorig) * c(ans0 / vsd)
inter <- drop(mean(Yorig) - sd(Yorig) * (vm / vsd) %*% ans0)
out\$intercept <- drop(mean(Yorig) - sd(Yorig) * (vm / vsd) %*% ans0)
out\$betahat <- sd(Yorig) * c(ans0 / vsd)
out\$best <- best
out\$Yinit <- Yinit

out\$yhat <- Xorig %*% sol + inter

out
}

############################################################
#
#  FUNCTIONS FOR CALCULATING NON-EUCLIDEAN MEANS AND DISTANCES
#  OF COVARIANCE MATRICES
#
############################################################
# Log Euclidean mean: Sigma_L

#==================================================================================
estLogEuclid <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
eS <- eigen(S[, , j], symmetric = TRUE)
sum <- sum +
weights[j] * eS\$vectors %*% diag(log(eS\$values)) %*% t(eS\$vectors) /
sum(weights)
}
ans <- sum
eL <- eigen(ans, symmetric = TRUE)
eL\$vectors %*% diag(exp(eL\$values)) %*% t(eL\$vectors)
}

#==================================================================================
estPowerEuclid <- function(S, weights = 1, alpha = 0.5) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
eS <- eigen(S[, , j], symmetric = TRUE)
sum <- sum +
weights[j] * eS\$vectors %*% diag(abs(eS\$values) ** alpha) %*% t(eS\$vectors) /
sum(weights)
}
ans <- sum
eL <- eigen(ans, symmetric = TRUE)
eL\$vectors %*% diag(abs(eL\$values) ** (1 / alpha)) %*% t(eL\$vectors)
}

# Riemannian (weighted mean) : Sigma_R
#==================================================================================
estLogRiem2 <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
check <- 9
tau <- 1
Hold <- 99999
mu <- estLogEuclid(S, weights)
while (check > 0.0000000001) {
ev <- eigen(mu, symmetric = TRUE)
logmu <- ev\$vectors %*% diag(log(ev\$values)) %*% t(ev\$vectors)
Hnew <- Re(Hessian2(S, mu, weights))
logmunew <- logmu + tau * Hnew
ev <- eigen(logmunew, symmetric = TRUE)
mu <- ev\$vectors %*% diag(exp(ev\$values)) %*% t(ev\$vectors)
check <- Re(Enorm(Hold) - Enorm(Hnew))
if (check < 0) {
tau <- tau / 2
check <- 999999
}
Hold <- Hnew
}
mu
}

# Hessian used in calculating Sigma_R
#==================================================================================
Hessian2 <- function(S, Sigma, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
ev0 <- eigen(Sigma, symmetric = TRUE)
shalf <-
ev0\$vectors %*% (diag(sqrt((ev0\$values)))) %*% t(ev0\$vectors)
sumit <- matrix(0, k, k)
for (i in 1:(M)) {
ev2 <- eigen(shalf %*% solve(S[, , i]) %*% shalf, symmetric = TRUE)
sumit <-
sumit + weights[i] * ev2\$vectors %*% diag (log((ev2\$values))) %*% t(ev2\$vectors) /
sum(weights)
}
- sumit
}

# Euclidean : Sigma_E
#==================================================================================
estEuclid <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
sum <- sum + S[, , j] * weights[j] / sum(weights)
}
sum
}

# Cholesky mean : Sigma_C
#==================================================================================
estCholesky <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
sum <- sum + t(chol(S[, , j])) * weights[j] / sum(weights)
}
cc <- sum
cc %*% t(cc)
}

#==================================================================================
ild_estSS <- function(S, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
H <- defh(k)
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
Q <- array(0, c(k + 1, k, M))
for (j in 1:M) {
Q[, , j] <- t(H) %*% (rootmat(S[, , j]))
}
ans <- procWGPA(
Q,
fixcovmatrix = diag(k + 1),
scale = FALSE,
reflect = TRUE,
sampleweights = weights
)
H %*% ans\$mshape %*% t(H %*% ans\$mshape)
}

#==================================================================================
ild_estShape <- function(S, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
H <- defh(k)
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
Q <- array(0, c(k + 1, k, M))
for (j in 1:M) {
Q[, , j] <- t(H) %*% (rootmat(S[, , j]))
}
ans <- procWGPA(
Q,
fixcovmatrix = diag(k + 1),
scale = TRUE,
reflect = TRUE,
sampleweights = weights
)
H %*% ans\$mshape %*% t(H %*% ans\$mshape)
}

#==================================================================================
estRiemLe <- function(S, weights) {
M <- dim(S)[3]
k <- dim(S)[1]
if (M != 2)
print("Sorry - Calculation not implemented for M>2 yet")
if (M == 2) {
P1 <- S[, , 1]
P2 <- S[, , 2]
detP1 <- prod(eigen(P1)\$values)
detP2 <- prod(eigen(P2)\$values)
P1 <- P1 / (detP1) ^ (1 / k)
P2 <- P2 / (detP2) ^ (1 / k)
P1inv <- solve(P1)
P12sq <- P1inv %*% P2 %*% P2 %*% P1inv
tem <- eigen(P12sq, symmetric = TRUE)
A2 <- tem\$vectors %*% diag(log(tem\$values)) %*% t(tem\$vectors)
logPs2 <- weights[2] * A2
tem2 <- eigen(logPs2, symmetric = TRUE)
Ps2 <- tem2\$vectors %*% diag(exp(tem2\$values)) %*% t(tem2\$vectors)
P12s <- P1 %*% Ps2 %*% P1
tem3 <- eigen(P12s, symmetric = TRUE)
P12sA <-
tem3\$vectors %*% diag(sqrt(tem3\$values)) %*% t(tem3\$vectors)
Ptildes <- (detP1 * (detP2 / detP1) ^ weights[2]) ^ (1 / k) * P12sA
Ptildes
}
}

##########distances#################################

#==================================================================================
distRiemPennec <- function(P1, P2) {
eig <- eigen(P1, symmetric = TRUE)
P1half <- eig\$vectors %*% diag(sqrt(eig\$values)) %*% t(eig\$vectors)
P1halfinv <- solve(P1half)
AA <- P1halfinv %*% P2 %*% P1halfinv
tem <- eigen(AA, symmetric = TRUE)
A2 <- tem\$vectors %*% diag(log(tem\$values)) %*% t(tem\$vectors)
dd <- Enorm(A2)
dd
}

#==================================================================================
distLogEuclidean <- function(P1, P2) {
eig <- eigen(P1, symmetric = TRUE)
logP1 <- eig\$vectors %*% diag(log(eig\$values)) %*% t(eig\$vectors)
tem <- eigen(P2, symmetric = TRUE)
logP2 <- tem\$vectors %*% diag(log(tem\$values)) %*% t(tem\$vectors)
dd <- Enorm(logP1 - logP2)
dd
}

#==================================================================================
distRiemannianLe <-
function(P1, P2) {
dd <- distRiemPennec(P1 %*% t(P1), P2 %*% t(P2)) / 2
dd
}

#==================================================================================
ild_distProcrustesSizeShape <-
function (P1, P2)
{
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% rootmat(P1)
Q2 <- t(H) %*% rootmat(P2)
ans <- sqrt(
centroid.size(Q1) ^ 2 + centroid.size(Q2) ^ 2 - 2 *
centroid.size(Q1) * centroid.size(Q2) * cos(riemdist(Q1,
Q2, reflect = TRUE))
)
ans
}

#==================================================================================
ild_distProcrustesFull <- function(P1, P2) {
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% rootmat(P1)
Q2 <- t(H) %*% rootmat(P2)
ans <- riemdist(Q1, Q2, reflect = TRUE)
ans
}

#==================================================================================
distPowerEuclidean <- function(P1, P2, alpha = 1 / 2) {
if (alpha != 0) {
eS <- eigen(P1, symmetric = TRUE)
Q1 <- eS\$vectors %*% diag(abs(eS\$values) ^ alpha) %*% t(eS\$vectors)
eS <- eigen(P2, symmetric = TRUE)
Q2 <- eS\$vectors %*% diag(abs(eS\$values) ^ alpha) %*% t(eS\$vectors)
dd <- Enorm(Q1 - Q2) / abs(alpha)
}
if (alpha == 0) {
dd <- distLogEuclidean(P1, P2)
}
dd
}

#==================================================================================
ild_distCholesky <- function(P1, P2) {
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% t(chol(P1))
Q2 <- t(H) %*% t(chol(P2))
ans <- Enorm(Q1 - Q2)
ans
}

#==================================================================================
distEuclidean <- function(P1, P2) {
ans <- Enorm(P1 - P2)
ans
}

##################

#==================================================================================
distcov <- function(S1,
S2 ,
method = "Riemannian",
alpha = 1 / 2) {
if (method == "Procrustes") {
dd <- distProcrustesSizeShape(S1, S2)
}
if (method == "ProcrustesShape") {
dd <- distProcrustesFull(S1, S2)
}
if (method == "Riemannian") {
dd <- distRiemPennec(S1, S2)
}
if (method == "Cholesky") {
dd <- distCholesky(S1, S2)
}
if (method == "Power") {
dd <- distPowerEuclidean(S1, S2, alpha)
}
if (method == "Euclidean") {
dd <- distEuclidean(S1, S2)
}
if (method == "LogEuclidean") {
dd <- distLogEuclidean(S1, S2)
}
if (method == "RiemannianLe") {
dd <- distRiemannianLe(S1, S2)
}
dd
}

#==================================================================================
estcov <-
function (S,
method = "Riemannian",
weights = 1,
alpha = 1 / 2,
MDSk = 2)
{
out <- list(
mean = 0,
sd = 0,
pco = 0,
eig = 0,
dist = 0
)
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
if (method == "Procrustes") {
dd <- estSS(S, weights)
}
if (method == "ProcrustesShape") {
dd <- estShape(S, weights)
}
if (method == "Riemannian") {
dd <- estLogRiem2(S, weights)
}
if (method == "Cholesky") {
dd <- estCholesky(S, weights)
}
if (method == "Power") {
dd <- estPowerEuclid(S, weights, alpha)
}
if (method == "Euclidean") {
dd <- estEuclid(S, weights)
}
if (method == "LogEuclidean") {
dd <- estLogEuclid(S, weights)
}
if (method == "RiemannianLe") {
dd <- estRiemLe(S, weights)
}
out\$mean <- dd
sum <- 0
for (i in 1:M) {
sum <-
sum + weights[i] * distcov(S[, , i], dd, method = method) ^ 2 / sum(weights)
}
out\$sd <- sqrt(sum)
dist <- matrix(0, M, M)
for (i in 2:M) {
for (j in 1:(i - 1)) {
dist[i, j] <- distcov(S[, , i], S[, , j], method = method)
dist[j, i] <- dist[i, j]
}
}
out\$dist <- dist
if (M > MDSk) {
ans <-
cmdscale(
dist,
k = MDSk,
eig = TRUE,
x.ret = TRUE
)
out\$pco <- ans\$points
out\$eig <- ans\$eig
if (MDSk > 2) {
shapes3d(out\$pco[, 1:min(MDSk, 3)], axes3 = TRUE)
}
if (MDSk == 2) {
plot(out\$pco,
type = "n",
xlab = "MDS1",
ylab = "MDS2")
text(out\$pco[, 1], out\$pco[, 2], 1:length(out\$pco[, 1]))
}
}
out
}

rootmat <- function(P1) {
eS <- eigen(P1, symmetric = TRUE)
if (min(eS\$values) < -0.001) {
print("Not positive-semi definite")
}
else{
Q1 <- eS\$vectors %*% diag(sqrt(abs(eS\$values))) %*% t(eS\$vectors)
Q1
}
}

##########################

#==================================================================================
shapes.cva <- function(X ,
groups ,
scale = TRUE,
tangentcoords = "residual",
ncv = 2) {
g <- dim(table (groups))
ans <- procGPA(X , tangentcoords=tangentcoords, scale = scale)

if (scale == TRUE)
pp <- (ans\$k - 1) * ans\$m - (ans\$m * (ans\$m - 1) / 2) - 1
if (scale == FALSE)
pp <- (ans\$k - 1) * ans\$m - (ans\$m * (ans\$m - 1) / 2)

pracdim <- min(pp,  ans\$n - g)
out <- lda(ans\$scores[, 1:pracdim] , groups)
print((out))
cv <- predict(out, dimen = ncv)\$x
if (dim(cv)[2] == 1) {
cv <- cbind(cv, rnorm(dim(cv)[1]) / 1000)
}
if (ncv == 2) {
eqscplot(cv,
type = "n",
xlab = "CV1",
ylab = "CV2")
text(cv, labels = groups)
}
if (ncv == 3) {
shapes3d(cv, color = groups, axes3 = TRUE)
}
cv
}

#==================================================================================
groupstack <- function(A1,
A2,
A3 = 0,
A4 = 0,
A5 = 0,
A6 = 0,
A7 = 0,
A8 = 0) {
out <- list(x = 0, groups = "")
dat <- abind(A1, A2)
group <- c(rep(1, times = dim(A1)[3]), rep(2, times = dim(A2)[3]))
if (is.array(A3)) {
dat <- abind(dat, A3)
group <- c(group, rep(3, times = dim(A3)[3]))
if (is.array(A4)) {
dat <- abind(dat, A4)
group <- c(group, rep(4, times = dim(A4)[3]))
if (is.array(A5)) {
dat <- abind(dat, A5)
group <- c(group, rep(5, times = dim(A5)[3]))
if (is.array(A6)) {
dat <- abind(dat, A6)
group <- c(group, rep(6, times = dim(A6)[3]))
if (is.array(A7)) {
dat <- abind(dat, A7)
group <- c(group, rep(7, times = dim(A7)[3]))
if (is.array(A8)) {
dat <- abind(dat, A8)
group <- c(group, rep(8, times = dim(A8)[3]))
}
}
}
}
}
}
out\$x <- dat
out\$groups <- group
out
}

###########################

#==================================================================================
procdist <- function(x, y, type = "full", reflect = FALSE) {
if (type == "full") {
out <- sin(riemdist(x, y, reflect = reflect))
}
if (type == "partial") {
out <- sqrt(2) * sqrt(abs(1 - cos(riemdist(x, y, reflect = reflect))))
}
if (type == "Riemannian") {
out <- riemdist(x, y, reflect = reflect)
}
if (type == "sizeandshape") {
out <- ssriemdist(x, y, reflect = reflect)
}
out
}

#==================================================================================
transformations <- function(Xrotated, Xoriginal) {
# outputs the translations, rotations and
# scalings for ordinary Procrustes rotation
# of each individual in Xoriginal to the
# Procrustes rotated individuals in Xrotated
X1 <- Xrotated
X2 <- Xoriginal
n <- dim(X1)[3]
m <- dim(X1)[2]
translation <- matrix(0, m, n)
scale <- rep(0, times = n)
rotation <- array(0, c(m, m, n))
for (i in 1:n) {
translation[, i] <- -apply(X2[, , i] - X1[, , i], 2, mean)
ans <- procOPA(X1[, , i], X2[, , i])
scale[i] <- ans\$s
rotation[, , i] <- ans\$R
}
out <- list(translation = 0,
scale = 0,
rotation = 0)
out\$translation <- translation
out\$scale <- scale
out\$rotation <- rotation
out
}

#==================================================================================
iglogl <- function(x , lam, nlam) {
gamma <- abs(x[1])
alpha <- gamma / mean(1 / lam[1:nlam])
ll <-
-(gamma + 1) * sum(log(lam[1:nlam])) - alpha * sum (1 / lam[1:nlam]) +
nlam * gamma * log(alpha) - nlam * lgamma (gamma)
- ll
}

#==================================================================================
procWGPA <-
function(x,
fixcovmatrix = FALSE,
initial = "Identity",
maxiterations = 10,
scale = TRUE,
reflect = FALSE,
prior = "Exponential",
diagonal = TRUE,
sampleweights = "Equal") {
X <- x
priorargument <- prior
alpha <- "not estimated"
gamma <- "not estimated"
k <- dim(X)[1]
n <- dim(X)[3]
m <- dim(X)[2]

if (initial[1] == "Identity") {
Sigmak <- diag(k)
}
else{
if (initial[1] == "Rawdata") {
tol <- 0.0000000001
if (m == 2) {
Sigmak <- diag(diag(var(t(X[, 1, ]))) + diag(var(t(X[, 2, ])))) / 2 + tol
}
if (m == 3) {
Sigmak <-
diag(diag(var(t(X[, 1, ]))) + diag(var(t(X[, 2, ]))) + diag(var(t(X[, 3, ])))) /
3 + tol
}
}
else
{
Sigmak <- initial
}
}

mu <- procGPA(X, scale = scale)\$mshape

#cat("Iteration 1 \n")

if (fixcovmatrix[1] != FALSE) {
Sigmak <- fixcovmatrix
}

ans <-
procWGPA1(
X,
mu,
metric = Sigmak,
scale = scale,
reflect = reflect,
sampleweights = sampleweights
)

if ((maxiterations > 1) && (fixcovmatrix[1] == FALSE)) {
ans0 <- ans
dif <- 999999
it <- 1
while ((dif > 0.00001) && (it < maxiterations)) {
it <- it + 1
if (it == 2) {
cat("Differences in norm of Sigma estimates... \n ")
}

if (prior[1] == "Identity") {
prior <- diag(k)
}

if (prior[1] == "Inversegamma") {
lam <- eigen(ans\$Sigmak)\$values
nlam <- min(c(n * m - m - 3, k - 3))
mu <- mean(1 / lam[1:(nlam)])
alpha <- 1 / mu
out <- nlm(iglogl,
p = c(1) ,
lam = lam,
nlam = nlam)
#print(out)
gamma <- abs(out\$estimate[1])
alpha <- gamma / mean(1 / lam[1:nlam])
newmetric <-
n * m / (n * m + 2 * (1 + gamma)) * (ans\$Sigmak + (2 * alpha / (n * m)) *
diag(k))
#dif2<-999999
#while (dif2> 0.000001){
#old<-alpha
#lam <- eigen(newmetric)\$values
#out <- nlm( iglogl, p=c(1) ,lam=lam, nlam=nlam)
#gamma <- abs(out\$estimate[1])
#alpha<- gamma/ mean(1/lam[1:nlam])
#newmetric <- n*m/(n*m+2*(1+gamma))*( ans\$Sigmak + (2*alpha/(n*m))*diag(k) )
#dif2<- abs(alpha- old)
#print(dif2)
#}
}

if (prior[1] == "Exponential") {
lam <- eigen(ans\$Sigmak)\$values
nlam <- min(c(n * m - m - 2, k - 2))
mu <- mean(1 / lam[1:(nlam)])
alpha <- 1 / mu
gamma <- 1
newmetric <-
n * m / (n * m + 2 * (2)) * (ans\$Sigmak + (2 * alpha / (n * m)) * diag(k))
#dif2<-999999
#while (dif2> 0.000001){
#old<-alpha
#newmetric <- n*m/(n*m+2*(2))*( ans\$Sigmak + (2*alpha/(n*m))*diag(k) )
#lam <- eigen(newmetric)\$values
#mu <- mean(1/lam[1:( nlam)])
#alpha <- 1/mu
#newmetric <- n*m/(n*m+2*(2))*( ans\$Sigmak + (2*alpha/(n*m))*diag(k) )
#dif2<- abs(alpha- old)
#}
}

if (is.double(prior[1])) {
newmetric <- (ans\$Sigmak + prior)
}

if (diagonal == TRUE) {
newmetric <- diag(diag(newmetric))
}

if (fixcovmatrix[1] != FALSE) {
newmetric <- fixcovmatrix
}

ans2 <-
procWGPA1(
X,
ans\$mshape,
metric = newmetric ,
scale = scale,
sampleweights = sampleweights
)
plotshapes(ans2\$rotated)
dif <- Enorm((ans\$Sigmak - ans2\$Sigmak))
ans <- ans2
cat(c(it, " ", dif, " \n"))
}
}
if ((priorargument[1] == "Exponential") ||
(priorargument[1] == "Inversegamma")) {
ans\$alpha <- alpha
ans\$gamma <- gamma
}
cat(" \n")
ans
}

#==================================================================================
procWGPA1 <- function(X,
mu,
metric = "Identity",
scale = TRUE,
reflect = FALSE,
sampleweights = "Equal") {
k <- dim(X)[1]
n <- dim(X)[3]
m <- dim(X)[2]

sum <- 0
for (i in 1:n) {
sum <- sum + centroid.size(X[, , i]) ** 2
}
size1 <- sqrt(sum)

if (sampleweights[1] == "Equal") {
sampleweights <- rep(1 / n, times = n)
}

if (length(sampleweights) != n) {
cat("Sample weight vector not of correct length \n")
}

if (metric[1] == "Identity") {
Sigmak <- diag(k)
}
else{
Sigmak <- metric
}

eig <- eigen(Sigmak, symmetric = TRUE)

Sighalf <-
eig\$vectors %*% diag (sqrt(abs(eig\$values))) %*% t(eig\$vectors)
Siginvhalf <-
eig\$vectors %*% diag(1 / sqrt(abs(eig\$values))) %*% t(eig\$vectors)
Siginv <- eig\$vectors %*% diag (1 / (eig\$values)) %*% t(eig\$vectors)

one <- matrix(rep(1, times = k), k, 1)

Xstar <- X

for (i in 1:n) {
Xstar[, , i] <- Xstar[, , i] - one %*% t(one) %*% Siginv %*% Xstar[, , i] /
c(t(one) %*% Siginv %*% one)
Xstar[, , i] <- Siginvhalf %*% Xstar[, , i]
}

mu <- mu - one %*% t(one) %*% Siginv %*% mu / c(t(one) %*% Siginv %*% one)

ans <- procGPA(Xstar, eigen2d = FALSE)
ans2 <- ans

dif3 <- 99999999
while (dif3 > 0.00001) {
for (i in 1:n) {
old <- mu
tem <-
procOPA(Siginvhalf %*% mu ,
Xstar[, , i],
scale = scale,
reflect = reflect)
Gammai <- tem\$R
betai <- tem\$s
#ci <- t(one)%*% Siginvhalf %*% X[,,i] %*% Gammai*betai/k
#Yi <- Sighalf%*% ans\$rotated[,,i] + Sighalf%*%one%*% ci
#Zi <- Yi - one %*% t(one)%*% Siginv %*% Yi / c( t(one)%*%Siginv%*%one )

Zi <- Sighalf %*% Xstar[, , i] %*% Gammai * betai
ans2\$rotated[, , i] <- Zi
}

sum2 <- 0
for (i in 1:n) {
sum2 <- sum2 + centroid.size(ans2\$rotated[, , i]) ** 2
}
size2 <- sqrt(sum2)

tem <- ans2\$mshape * 0
for (i in 1:n) {
ans2\$rotated[, , i] <- ans2\$rotated[, , i] * size1 / size2
tem <- tem + ans2\$rotated[, , i] * sampleweights[i] / sum(sampleweights)
}

mu <- tem
dif3 <- riemdist(old,  mu)
}

z <- ans2
z\$mshape <- tem

tan <- z\$rotated[, 1,] - z\$mshape[, 1]
for (i in 2:m) {
tan <- rbind(tan, z\$rotated[, i,] - z\$mshape[, i])
}
pca <- prcomp1(t(tan))
z\$tan <- tan
npc <- 0
for (i in 1:length(pca\$sdev)) {
if (pca\$sdev[i] > 1e-07) {
npc <- npc + 1
}
}
z\$scores <- pca\$x
z\$rawscores <- pca\$x
z\$stdscores <- pca\$x
for (i in 1:npc) {
z\$stdscores[, i] <- pca\$x[, i] / pca\$sdev[i]
}
z\$pcar <- pca\$rotation
z\$pcasd <- pca\$sdev
z\$percent <- z\$pcasd ^ 2 / sum(z\$pcasd ^ 2) * 100

size <- rep(0, times = n)
rho <- rep(0, times = n)
x <- X
size <- apply(x, 3, centroid.size)
rho <- apply(x, 3, y <- function(x) {
riemdist(x, z\$mshape)
})

z\$rho <- rho
z\$size <- size
z\$rmsrho <- sqrt(mean(rho ^ 2))
z\$rmsd1 <- sqrt(mean(sin(rho) ^ 2))

z\$k <- k
z\$m <- m
z\$n <- n

tem <- matrix(0, k, k)

for (i in 1:n) {
tem <-
tem +  (z\$rotated[, , i] - z\$mshape) %*% t((z\$rotated[, , i] - z\$mshape))
}
tem <- tem / (n * m)
z\$Sigmak <-  tem

return(z)
}

#==================================================================================
testshapes <-
function(A,
B,
resamples = 1000,
replace = TRUE,
scale = TRUE) {
if (replace == TRUE) {
out <- bootstraptest(A, B, resamples = resamples, scale = scale)
}
if (replace == FALSE) {
out <- permutationtest(A, B, nperms = resamples, scale = scale)
}
out
}

#==================================================================================
testmeanshapes <-
function(A,
B,
resamples = 1000,
replace = FALSE,
scale = TRUE) {
if (replace == TRUE) {
out <- bootstraptest(A, B, resamples = resamples, scale = scale)
}
if (replace == FALSE) {
out <- permutationtest(A, B, nperms = resamples, scale = scale)
}

if (resamples > 0) {
aa <- list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0
)

aa\$H <- out\$H
aa\$H.pvalue <- out\$H.pvalue
aa\$H.table.pvalue <- out\$H.table.pvalue
aa\$G <- out\$G
aa\$G.pvalue <- out\$G.pvalue
aa\$G.table.pvalue <- out\$G.table.pvalue
aa\$J <- out\$J
aa\$J.pvalue <- out\$J.pvalue
aa\$J.table.pvalue <- out\$J.table.pvalue

}

if (resamples == 0) {
aa <- list(
H = 0,
H.table.pvalue = 0,
G = 0,
G.table.pvalue = 0,
J = 0,
J.table.pvalue = 0
)

aa\$H <- out\$H
aa\$H.table.pvalue <- out\$H.table.pvalue
aa\$G <- out\$G
aa\$G.table.pvalue <- out\$G.table.pvalue
aa\$J <- out\$J
aa\$J.table.pvalue <- out\$J.table.pvalue

}

aa
}

#==================================================================================
permutationtest2 <- function (A, B, nperms = 1000, scale = scale)
{
A1 <- A
A2 <- B
mdim <- dim(A1)[2]
B <- nperms
nsam1 <- dim(A1)[3]
nsam2 <- dim(A2)[3]
pool <-
procGPA(
abind (A1, A2) ,
scale = scale,
tangentcoords = "partial",
pcaoutput = FALSE
)

tempool <- pool
for (i in 1:(nsam1 + nsam2)) {
tempool\$tan[, i] <- pool\$tan[, i] / Enorm(pool\$tan[, i]) * pool\$rho[i]
}
pool <- tempool

permpool <- pool
Gtem <- Goodall(pool, nsam1, nsam2)
Htem <- Hotelling(pool, nsam1, nsam2)
Jtem <- James(pool, nsam1, nsam2, table = TRUE)
Ltem <- Lambdamin(pool, nsam1, nsam2)
Gumc <- Gtem\$F
Humc <- Htem\$F
Jumc <- Jtem\$Tsq
Lumc <- Ltem\$lambdamin
Gtabpval <- Gtem\$pval
Htabpval <- Htem\$pval
Jtabpval <- Jtem\$pval
Ltabpval <- Ltem\$pval
if (B > 0) {
Apool <- array(0, c(dim(A1)[1], dim(A1)[2], dim(A1)[3] +
dim(A2)[3]))
Apool[, , 1:nsam1] <- A1
Apool[, , (nsam1 + 1):(nsam1 + nsam2)] <- A2
out <-
list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0
)
Gu <- rep(0, times = B)
Hu <- rep(0, times = B)
Ju <- rep(0, times = B)
Lu <- rep(0, times = B)
cat("Permutations - sampling without replacement: ")
cat(c("No of permutations = ", B, "\n"))
for (i in 1:B) {
if (i / 100 == trunc(i / 100)) {
cat(c(i, " "))
}
select <- sample(1:(nsam1 + nsam2))
permpool\$tan <- pool\$tan[, select]
Gu[i] <- Goodall(permpool, nsam1, nsam2)\$F
Hu[i] <- Hotelling(permpool, nsam1, nsam2)\$F
Ju[i] <- James(permpool, nsam1, nsam2)\$Tsq
Lu[i] <-
Lambdamin(permpool, nsam1, nsam2)\$lambdamin
}
Gu <- sort(Gu)
numbig <- length(Gu[Gumc < Gu])
pvalG <- (1 + numbig) / (B + 1)
Lu <- sort(Lu)
numbig <- length(Lu[Lumc < Lu])
pvalL <- (1 + numbig) / (B + 1)
Ju <- sort(Ju)
numbig <- length(Ju[Jumc < Ju])
pvalJ <- (1 + numbig) / (B + 1)
Hu <- sort(Hu)
numbig <- length(Hu[Humc < Hu])
pvalH <- (1 + numbig) / (B + 1)
cat(" \n")
out\$Hu <- Hu
out\$Ju <- Ju
out\$Gu <- Gu
out\$Lu <- Lu
out\$H <- Humc
out\$H.pvalue <- pvalH
out\$H.table.pvalue <- Htabpval
out\$J <- Jumc
out\$J.pvalue <- pvalJ
out\$J.table.pvalue <- Jtabpval
out\$G <- Gumc
out\$G.pvalue <- pvalG
out\$G.table.pvalue <- Gtabpval
out\$L <- Lumc
out\$L.pvalue <- pvalL
out\$L.table.pvalue <- Ltabpval
}
if (B == 0) {
out <- list(
H = 0,
H.table.pvalue = 0,
G = 0,
G.table.pvalue = 0
)
out\$H <- Humc
out\$H.table.pvalue <- Htabpval
out\$J <- Jumc
out\$J.table.pvalue <- Jtabpval
out\$G <- Gumc
out\$G.table.pvalue <- Gtabpval
out\$L <- Lumc
out\$L.table.pvalue <- Ltabpval
}
out
}

#==================================================================================
bootstraptest <- function (A,
B,
resamples = 200,
scale = TRUE)
{
A1 <- A
A2 <- B
mdim <- dim(A1)[2]
B <- resamples
nsam1 <- dim(A1)[3]
nsam2 <- dim(A2)[3]
pool <-
procGPA(
abind (A1, A2) ,
scale = scale ,
tangentcoords = "partial",
pcaoutput = FALSE
)
tempool <- pool
for (i in 1:(nsam1 + nsam2)) {
tempool\$tan[, i] <- pool\$tan[, i] / Enorm(pool\$tan[, i]) * pool\$rho[i]
}
pool <- tempool
bootpool <- pool
Gtem <- Goodall(pool, nsam1, nsam2)
Htem <- Hotelling(pool, nsam1, nsam2)
Jtem <- James(pool, nsam1, nsam2, table = TRUE)
Ltem <- Lambdamin(pool, nsam1, nsam2)
Gumc <- Gtem\$F
Humc <- Htem\$F
Jumc <- Jtem\$Tsq
Lumc <- Ltem\$lambdamin
Gtabpval <- Gtem\$pval
Htabpval <- Htem\$pval
Jtabpval <- Jtem\$pval
Ltabpval <- Ltem\$pval
if (B > 0) {
Apool <- array(0, c(dim(A1)[1], dim(A1)[2], dim(A1)[3] +
dim(A2)[3]))
Apool[, , 1:nsam1] <- A1
Apool[, , (nsam1 + 1):(nsam1 + nsam2)] <- A2
out <-
list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0
)
Gu <- rep(0, times = B)
Hu <- rep(0, times = B)
Ju <- rep(0, times = B)
Lu <- rep(0, times = B)
pool2 <- pool
pool2\$tan[, 1:nsam1] <-
pool\$tan[, 1:nsam1] - apply(pool\$tan[, 1:nsam1], 1, mean)
pool2\$```