R/shapes.R

Defines functions dif.old dis del cnt3 close1 bgpa add MDSshape objfun4 objfun frechet permutationtest oneFone ild_kendall.shpv loneFone isologdens isodens loglikeiso2 loglikeiso objfuniso isomle banner4 banner1 plot3Dpca plot3Dmean plot3Ddata.static plot3Ddata tanpreshape ild_tanfigurefull ild_tanfigure st sigmacov rotateaxes riemdist.mD riemdist.complex riemdist ssriemdist relwarps reassqpr realtocomplex read.in read.array project procGPA testmeanshapes.old procrustes2d procdistreflect prinwscoregrids ild_preshapetoicon ild_preshape.mat ild_preshape.mD ild_preshape plotrelwarp plotproc plotprinwarp plotpairscores pointsPDMnoaxis3 plotPDMnoaxis plotPDMbook plotPDM3 plotPDM2 plotPDM plot2rwscores partialwarps partialwarpgrids partial.procdist ild_Enorm movie makearray mahpreshapedist linegrid isotropy.test genpower full.procdist ild_defh complextoreal ild_centroid.size.mD ild_centroid.size.complex ild_centroid.size cbevectors cbevec bookstein.shpv.complex bookstein.shpv bookstein2d shaperw Vmat Vinv V I2mat Goodalltest Goodall2D BoxM plotshapes defplotsize2 procOPA defplotsize3 resampletest MGM shapes3d abind as.3d rigidbody rotatexyz James Hotelling Goodall testmeanshapes testshapes procWGPA1 procWGPA iglogl transformations procdist groupstack shapes.cva rootmat distcov distEuclidean ild_distCholesky distPowerEuclidean ild_distProcrustesFull distRiemannianLe distLogEuclidean distRiemPennec estRiemLe ild_estShape ild_estSS estCholesky estEuclid Hessian2 estLogRiem2 estPowerEuclid estLogEuclid pedreg ped project.subsphere col2RGB centroid.size defh multiply_by_helmert_implicitly_3d uji3_centroid.size uji2_centroid.size multiply_by_transpose_of_helmert_implicitly multiply_by_transpose_of_helmert_explicitly multiply_by_transpose_of_helmert multiply_by_helmert_implicitly multiply_by_helmert_explicitly multiply_by_helmert uji_kendall.shpv uji_tanfigurefull uji_tanfigure uji_preshape.mat uji_preshape.mD uji_centroid.size.mD uji_centroid.size.complex uji_estShape uji_estSS uji_distCholesky uji_distProcrustesSizeShape uji_distProcrustesFull uji_Enorm uji_defh uji_centroid.size uji_preshape fort.ROTATION vec1 sim1 sh sgpa rgpa msh graf ftrsq fos.REFLECT fort.ROTATEANDREFLECT fopa fgpa.rot fgpa fgpa.singleiteration fcnt fcel fJ Plot3D plotshapes3d.pns shape.pcscores.partial shape.pcscores tangent.coords.partial Procrustes.dist.full rot.mat pc2sphere pcscore2sphere sphere2pcscore Enormalize tr sphereFit sphere.jac sphere.res sphere.obj flipud0 repmat mod get.data.subsphere get.prinarc.subsphere get.prinarc get.prinarc.value trans.subsphere UNIFORMdirections randvonMisesFisherm PNSs2e PNSe2s geodmeanS1 vMFtest LRTpval getSubSphere objfn LogNPd ExpNPd rotMat pns.pc pns4pc pns plot3darcs preshape2shape tangentcoords.partial.inv

Documented in abind add as.3d banner1 banner4 bgpa bookstein2d bookstein.shpv bookstein.shpv.complex BoxM cbevec cbevectors centroid.size close1 cnt3 complextoreal defh defplotsize2 defplotsize3 del dif.old dis distcov distEuclidean distLogEuclidean distPowerEuclidean distRiemannianLe distRiemPennec estCholesky estEuclid estLogEuclid estLogRiem2 estPowerEuclid estRiemLe fcel fcnt fgpa fgpa.rot fgpa.singleiteration fJ fopa fort.ROTATEANDREFLECT fort.ROTATION fos.REFLECT frechet ftrsq full.procdist genpower Goodall Goodall2D Goodalltest graf groupstack Hessian2 Hotelling I2mat iglogl ild_centroid.size ild_centroid.size.complex ild_centroid.size.mD ild_defh ild_distCholesky ild_distProcrustesFull ild_Enorm ild_estShape ild_estSS ild_kendall.shpv ild_preshape ild_preshape.mat ild_preshape.mD ild_preshapetoicon ild_tanfigure ild_tanfigurefull isodens isologdens isomle isotropy.test James linegrid loglikeiso loglikeiso2 loneFone mahpreshapedist makearray MDSshape MGM movie msh multiply_by_helmert multiply_by_helmert_explicitly multiply_by_helmert_implicitly multiply_by_helmert_implicitly_3d multiply_by_transpose_of_helmert multiply_by_transpose_of_helmert_explicitly multiply_by_transpose_of_helmert_implicitly objfun objfun4 objfuniso oneFone partial.procdist partialwarpgrids partialwarps ped pedreg permutationtest plot2rwscores plot3darcs plot3Ddata plot3Ddata.static plot3Dmean plot3Dpca plotpairscores plotPDM plotPDM2 plotPDM3 plotPDMbook plotPDMnoaxis plotprinwarp plotproc plotrelwarp plotshapes pns pns4pc pointsPDMnoaxis3 preshape2shape prinwscoregrids procdist procdistreflect procGPA procOPA procrustes2d procWGPA procWGPA1 project read.array read.in realtocomplex reassqpr relwarps resampletest rgpa riemdist riemdist.complex riemdist.mD rigidbody rootmat rotateaxes rotatexyz sgpa sh shaperw shapes3d shapes.cva sigmacov sim1 ssriemdist st tangentcoords.partial.inv tanpreshape testmeanshapes testmeanshapes.old testshapes transformations uji2_centroid.size uji3_centroid.size uji_centroid.size uji_centroid.size.complex uji_centroid.size.mD uji_defh uji_distCholesky uji_distProcrustesFull uji_distProcrustesSizeShape uji_Enorm uji_estShape uji_estSS uji_kendall.shpv uji_preshape uji_preshape.mat uji_preshape.mD uji_tanfigure uji_tanfigurefull V vec1 Vinv Vmat

#-----------------------------------------------------------------------
#
# Statistical shape analysis routines
# written by Ian Dryden in R  (see http://cran.r-project.org)
# (c) Ian Dryden
#     University of Nottingham. version 1.2.6
#                    2003-2021
#
# 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))
  # PNS mean
  PNSmean = pns.out$PNS$mean
  GPAout = pns.out$GPAout
  
  # c * sd : lower and upper along principal arcs. lower -> 0 -> upper
  # total 2 * nn + 1 configurations
  #    : nn configurations from negative to 0
  #      1 configuration at 0
  #      nn configurations from 0 to positive
  #default c = 1
  {
    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))
  # spherical data to PC scores
  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)
  # PC scores to tangent coordinates
  U1 = matrix(0, npts, ncol(GPAout$pcar))
  U2 = matrix(0, npts, ncol(GPAout$pcar))
  U3 = matrix(0, npts, ncol(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, ncol(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), ncol(GPAout$pcar))
  tan.lu.arc2 = matrix(0, nrow(lu.arc2), ncol(GPAout$pcar))
  tan.lu.arc3 = matrix(0, nrow(lu.arc3), ncol(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]
    }
  }
  # tangent coordinates to preshapes
  # preshapes to shapes
  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)
  {
    shapes.arc1[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U1[i, ], p = H %*% GPAout$mshape, R = diag(m)))
    shapes.arc2[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U2[i, ], p = H %*% GPAout$mshape, R = diag(m)))
    shapes.arc3[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U3[i, ], p = H %*% GPAout$mshape, R = diag(m)))
  }
  shapes.PNSmean = preshape2shape(tangentcoords.partial.inv(v = U.mean, 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))
  {
    shapes.lu.arc1[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc1[i, ], p = H %*% GPAout$mshape, R = diag(m)))
    shapes.lu.arc2[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc2[i, ], p = H %*% GPAout$mshape, R = diag(m)))
    shapes.lu.arc3[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc3[i, ], p = H %*% GPAout$mshape, R = diag(m)))
  }
  
  
  #convert PCs to configuration space if not already in place 
  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
  }
  
  
  #------------------------#
  #                        #
  # PLOT TYPE I            #
  #    : full plot         #
  #                        #
  #------------------------#
 
  
   
   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)
     
    # GPA mean + PC
    plot3d(GPAout$mshape, type = "s", col = rainbow(k), size = 1, add = TRUE)
    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])   
                                                       
    spheres3d( pcu, radius = 0.004, color = "black")
    spheres3d( pcl, radius = 0.004, color = "grey")
    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)
 
    # principal arcs in landmark space
 
    # PNS mean in landmark space
    plot3d(shapes.PNSmean, type = "s", col = rainbow(k), size = 1.3, add = TRUE)
    lines3d(shapes.PNSmean,lwd=5,col=rainbow(k))
  
    
    # lower and upper
    for (i in 1:k)
    {
      lines3d(t(shapes.lu.arc[i, , ]), col = rainbow(k)[i], lwd = 1,lty=2)
      # farthest negative from 0
      spheres3d(head(t(shapes.lu.arc[i, , ]), 1), radius = 0.004, color = "black")
      if (i>1){
        lines3d( (shapes.lu.arc[(i-1):i,,1]) ,col="black" )
      }
      
      #   # farthest positive from 0
      spheres3d(tail(t(shapes.lu.arc[i, , ]), 1), radius = 0.004, 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 = 0, n.pc = 3) 
{
  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)
  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)
  print("Radii of spheres")
  print(pns.out$PNS$radii)
  pns.out$percent = pns.out$percent * sum(GPAout$percent[1:n.pc])/100
  pns.out$GPAout = GPAout
  pns.out$spheredata = spheredata
  return(pns.out)
}

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

  GPAout = procGPA(x = x, scale = TRUE, reflect = FALSE, tol1=1e-5,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))
}







#==================================================================================
# 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 = 0)
{
   n = ncol(x)
   k = nrow(x)
   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
   cat("Message from pns() : dataset is on ", d, "-sphere. \n", sep = "")
   if (nullspdim > 0)
   {
      cat(" .. found null space of dimension ",
          nullspdim,
          ", to be trivially reduced. \n",
          sep = "")
   }
   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)
         cat(d - i + 1,
             "-sphere to ",
             d - i,
             "-sphere, by ",
             "NULL space \n",
             sep = "")
      }
   }
   if (sphere.type == "seq.test")
   {
      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
               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
                  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
                  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
            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)
         {
            plot3d(
               x = t(currentSphere),
               xlab = "",
               ylab = "",
               zlab = "",
               xlim = c(-1, 1),
               ylim = c(-1, 1),
               zlim = c(-1, 1),
               box = FALSE,
               axes = FALSE,
               aspect = "iso"
            )
            axis3d(
               edge = 'x',
               labels = TRUE,
               tick = TRUE,
               at = c(-1, 1),
               pos = c(NA, 0, 0)
            )
            axis3d(
               edge = 'y',
               labels = TRUE,
               tick = TRUE,
               at = c(-1, 1),
               pos = c(0, NA, 0)
            )
            axis3d(
               edge = 'z',
               labels = TRUE,
               tick = TRUE,
               at = c(-1, 1),
               pos = c(0, 0, NA)
            )
         }
         
         if (nrow(currentSphere) == 2)
         {
            points3d(t(cur.proj), col = "red", size = 2)
         }
         
      }
   } else if (sphere.type == "BIC") {
      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)
         cat("BICsm: ", BICsmall, ", BICgr: ", BICgreat, "\n", sep = "")
         if (BICsmall > BICgreat)
         {
            center = center.g
            r = r.g
            cat(d - i + 1,
                "-sphere to ",
                d - i,
                "-sphere, by ",
                "GREAT sphere, BIC \n",
                sep = "")
         } else {
            center = center.s
            r = r.s
            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
         NestedSphere = rotMat(center) %*% currentSphere
         currentSphere = NestedSphere[1:(k - i),] /
            repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
                                  2), nrow = 1), k - i, 1)
      }
   } 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
         NestedSphere = rotMat(center) %*% currentSphere
         currentSphere = NestedSphere[1:(k - i),] /
            repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
                                  2), nrow = 1), k - i, 1)
      }
   } 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
            NestedSphere = rotMat(center) %*% currentSphere
            currentSphere = NestedSphere[1:(k - i),] /
               repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
                                     2), nrow = 1), k - i, 1)
         }
      }
      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
            NestedSphere = rotMat(center) %*% currentSphere
            currentSphere = NestedSphere[1:(k - i),] /
               repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
                                     2), nrow = 1), k - i, 1)
         }
      }
   } else {
      print("!!! Error from pns():")
      print("!!! sphere.type must be 'seq.test', 'small', 'great', 'BIC', or 'bi.sphere'")
      print("!!!   Terminating execution     ")
      return(NULL)
   }
   S1toRadian = atan2(currentSphere[2,], currentSphere[1,])
   meantheta = geodmeanS1(S1toRadian)$geodmean
   orthaxis[[d]] = meantheta
   resmat[d,] = mod(S1toRadian - meantheta + pi, 2 * pi) - pi
   
   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(
      cos(S1toRadian[which.max(l)]),
      sin(S1toRadian[which.max(l)]),
      pch = 4,
      cex = 3,
      col = "blue"
   )
   points(
      cos(S1toRadian[which.min(l)]),
      sin(S1toRadian[which.min(l)]),
      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")
   }
   radii = 1
   for (i in 1:(d - 1))
   {
      radii = c(radii, prod(sin(dist[1:i])))
   }
   resmat = flipud0(repmat(matrix(radii, ncol = 1), 1, n) * resmat)
   
   PNS = list()
   PNS$radii = radii
   PNS$orthaxis = orthaxis
   PNS$dist = dist
   PNS$pvalues = pvalues
   PNS$ratio = ratio
   PNS$basisu = NULL
   PNS$mean = c(PNSe2s(matrix(0, d, 1), PNS))
   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"
   }
   varPNS = apply(abs(resmat) ^ 2, 1, sum) / n
   total = sum(varPNS)
   propPNS = varPNS / total * 100
   return(list(
      resmat = resmat,
      PNS = PNS,
      percent = propPNS
   ))
}

#==================================================================================
pns4pc = function(x,
                  sphere.type = "seq.test",
                  alpha = 0.1,
                  R = 100,
                  nlast.small.sphere = 0,
                  n.pc = 2)
{
   if (n.pc < 2)
   {
      stop("Error: n.pc should be >= 2.")
   }
   out = pc2sphere(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
   U = princomp(t(x))$loadings[,]
   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
   radialdistances = acos(t(center.s) %*% x)
   xi_sample = mean(radialdistances) / sd(radialdistances)
   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
      radialdistances = acos(t(center.s) %*% rdata)
      xi_vec[r] = mean(radialdistances) / sd(radialdistances)
   }
   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,])),
            cos(NSradius[1] + res[2,]))
   if (dm > 2)
   {
      for (i in 1:(dm - 2))
      {
         T = t(rotMat(NSOrthaxis[[i + 1]])) %*%
            rbind(repmat(sin(NSradius[i + 1] + matrix(
               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)
   }
   S1toRadian = atan2(currentSphere[2,], currentSphere[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 = 300)
   )
   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)
}

#==================================================================================
pc2sphere = function(x, n.pc)
{
   k = dim(x)[1]
   m = dim(x)[2]
   n = dim(x)[3]
   if (n < ((k - 1) * m))
   {
      stop("Error: n 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,
         add = TRUE
      )
   }
   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,
      add = TRUE
   )
   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)
            gradL <- gradL + matrix(
               lambda / 2 * pars / Enorm(pars) * sum(abs(pars)) / GM +
                  lambda / 2 * sign(pars) * Enorm(pars) / GM ,
               p,
               1
            )
            
            
            c(gradL)
         }
         
         
         
         
         if (NUM_METHOD == 1) {
            repeat {
               #,ndeps=1e-3,factr=1e-5,pgtol=1e-5
               res2 <-
                  optim(
                     par = x0,
                     fn = pedfun,
                     gr = pedgrad,
                     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,
                     gr = pedgrad,
                     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,
               add = 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,
                       ncv = 2) {
   g <- dim(table (groups))
   ans <- procGPA(X , 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 = 3)$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$tan[, (nsam1 + 1):(nsam1 + nsam2)] <-
         pool$tan[, (nsam1 + 1):(nsam1 + nsam2)] -
         apply(pool$tan[, (nsam1 + 1):(nsam1 + nsam2)], 1, mean)
      
      
      cat("Bootstrap - sampling with replacement within each group under H0: ")
      cat(c("No of resamples = ", B, "\n"))
      for (i in 1:B) {
         if (i / 100 == trunc(i / 100)) {
            cat(c(i, " "))
         }
         select1 <- sample(1:nsam1, replace = TRUE)
         select2 <- sample((nsam1 + 1):(nsam1 + nsam2), replace = TRUE)
         bootpool$tan <- pool2$tan[, c(select1, select2)]
         Gu[i] <- Goodall(bootpool, nsam1, nsam2)$F
         Hu[i] <- Hotelling(bootpool, nsam1, nsam2)$F
         Ju[i] <- James(bootpool, nsam1, nsam2)$Tsq
         Lu[i] <- Lambdamin(bootpool, nsam1, nsam2)$lambdamin
      }
      Gu <- sort(Gu)
      numbig <- length(Gu[Gumc < Gu])
      pvalG <- (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)
      numbig <- length(Lu[Lumc < Lu])
      pvalL <- (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,
            J = 0,
            J.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
}

#==================================================================================
Lambdamin <-
   function (pool, n1, n2, p = 0)
   {
      censiz <- centroid.size(pool$mshape)
      tan1 <- pool$tan[, 1:n1]
      tan2 <- pool$tan[, (n1 + 1):(n1 + n2)]
      kt <- dim(tan1)[1]
      n <- n1 + n2
      k <- pool$k
      m <- pool$m
      if (p == 0) {
         p <- min(k * m - (m * (m - 1)) / 2 - 1 - m, n1 + n2 - 2)
      }
      
      
      HH <- diag(k)
      
      
      mu1 <- pool$mshape
      
      if (dim(tan1)[1] == k * m - m) {
         HH <- defh(k - 1)
         
         mu1 <- preshape(pool$mshape)
      }
      
      
      if (m == 2) {
         mu <- c(mu1[, 1], mu1[, 2])
         
      }
      if (m == 3) {
         mu <- c(mu1[, 1], mu1[, 2], mu1[,
                                         3])
         
      }
      
      
      
      
      
      dd <- kt
      X1 <- tan1 * 0
      X2 <- tan2 * 0
      S1 <- matrix(0, dd, dd)
      S2 <- matrix(0, dd, dd)
      for (i in 1:n1) {
         X1[, i] <- (mu + tan1[, i]) / Enorm(mu + tan1[, i])
         S1 <- S1 + X1[, i] %*% t(X1[, i])
      }
      for (i in 1:n2) {
         X2[, i] <- (mu + tan2[, i]) / Enorm(mu + tan2[, i])
         S2 <- S2 + X2[, i] %*% t(X2[, i])
      }
      
      sumx1 <- 0
      sumx2 <- 0
      for (i in 1:n1) {
         sumx1 <- sumx1 + X1[, i]
      }
      for (i in 1:n2) {
         sumx2 <- sumx2 + X2[, i]
      }
      
      
      sum1 <- apply(X1, 1, sum)
      sum2 <- apply(X2, 1, sum)
      mean1 <- sum1 / Enorm(sum1)
      mean2 <- sum2 / Enorm(sum2)
      
      
      bb1 <- mean1[1:(dd - 1)]
      cc1 <- mean1[dd]
      bb2 <- mean2[1:(dd - 1)]
      cc2 <- mean2[dd]
      A1 <- cc1 / abs(cc1) * diag(dd - 1) - cc1 / (abs(cc1) + cc1 ^ 2) *
         bb1 %*% t(bb1)
      M1 <- cbind(A1,-bb1)
      A1 <- cc2 / abs(cc2) * diag(dd - 1) - cc2 / (abs(cc2) + cc2 ^ 2) *
         bb2 %*% t(bb2)
      M2 <- cbind(A1,-bb2)
      G1 <- matrix(0, dd - 1, dd - 1)
      G2 <- matrix(0, dd - 1, dd - 1)
      for (iu in 1:(dd - 1)) {
         for (iv in iu:(dd - 1)) {
            G1[iu, iv] <- G1[iu, iv] + t((t(M1))[, iu]) %*% S1 %*%
               (t(M1))[, iv]
            G1[iv, iu] <- G1[iu, iv]
            G2[iu, iv] <- G2[iu, iv] + t((t(M2))[, iu]) %*% S2 %*%
               (t(M2))[, iv]
            G2[iv, iu] <- G2[iu, iv]
         }
      }
      G1 <- G1 / n1 / Enorm(sumx1 / n1) ^ 2
      G2 <- G2 / n2 / Enorm(sumx2 / n2) ^ 2
      #   eva1 <- eigen(G1, symmetric = TRUE, EISPACK = TRUE)
      eva1 <- eigen(G1, symmetric = TRUE)
      pcar1 <- eva1$vectors[, 1:p]
      pcasd1 <- sqrt(abs(eva1$values[1:p]))
      #    eva2 <- eigen(G2, symmetric = TRUE, EISPACK = TRUE)
      eva2 <- eigen(G2, symmetric = TRUE)
      pcar2 <- eva2$vectors[, 1:p]
      pcasd2 <- sqrt(abs(eva2$values[1:p]))
      if ((pcasd1[p] < 1e-06) || (pcasd2[p] < 1e-06)) {
         offset <- 1e-06
         cat("*")
         pcasd1 <- sqrt(pcasd1 ^ 2 + offset)
         pcasd2 <- sqrt(pcasd2 ^ 2 + offset)
      }
      Ahat1 <-
         n1 * t(M1) %*% (pcar1 %*% diag(1 / pcasd1 ^ 2) %*% t(pcar1)) %*%
         M1
      Ahat2 <-
         n2 * t(M2) %*% (pcar2 %*% diag(1 / pcasd2 ^ 2) %*% t(pcar2)) %*%
         M2
      
      
      Ahat <- (Ahat1 +  Ahat2)
      #    eva <- eigen(Ahat, symmetric = TRUE, EISPACK = TRUE)
      eva <- eigen(Ahat, symmetric = TRUE)
      
      
      
      lambdamin <- eva$values[p + 1]
      pval <- 1 - pchisq(lambdamin, p)
      #print(lambdamin)
      #print(pval)
      z <- list()
      z$pval <- pval
      z$df <- p
      z$lambdamin <- lambdamin
      return(z)
   }

#==================================================================================
Goodall <- function(pool , n1, n2, p = 0) {
   tan1 <- pool$tan[, 1:n1]
   tan2 <- pool$tan[, (n1 + 1):(n1 + n2)]
   kt <- dim(tan1)[1]
   n <- n1 + n2
   k <- pool$k
   m <- pool$m
   if (p == 0) {
      p <- min(k * m - (m * (m - 1)) / 2 - 1 - m, n1 + n2 - 2)
   }
   top <-  Enorm(apply(tan1, 1, mean) - apply(tan2, 1, mean)) ** 2
   bot <-
      sum(diag(var(t(tan1)))) * (n1 - 1) +  sum(diag(var(t(tan2)))) * (n2 - 1)
   Fstat <- ((n1 + n2 - 2) / (1 / n1 + 1 / n2) * top) / bot
   pval <- 1 - pf(Fstat, p, (n1 + n2 - 2) * p)
   z <- list()
   z$F <- Fstat
   z$pval <- pval
   z$df1 <- p
   z$df2 <- (n1 + n2 - 2) * p