New functions

.Rbuildignore

@@ -3,5 +3,3 @@
 ^\.github$
 ^LICENSE\.md$
 ^coverage_check\.R$
-^vignettes/annulus_demo\.Rmd$
-^vignettes/torus_demo\.Rmd$

.github/workflows/check-standard.yaml

@@ -39,7 +39,7 @@ jobs:
         with:
           r-version: ${{ matrix.config.r }}
 
-      - uses: r-lib/actions/setup-pandoc@v1
+      - uses: r-lib/actions/setup-pandoc@v2
 
       - name: Query dependencies
         run: |

DESCRIPTION

@@ -4,7 +4,8 @@ Version: 0.0.0.9000
 Authors@R: 
     person("Emily", "Winn-Nunez", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-6759-5406"))
-Description: This package allows users to generate alpha shapes from a distribution
+Description: Shape statistics is a growing field. This package allows users to 
+              generate alpha shapes from a distribution
               in two and three dimensions. It also allows users to take a data
               set of shapes and generate new alpha shapes from the distribution
               using non-black box methods.
@@ -17,7 +18,6 @@ Imports:
     stats,
     Rvcg,
     TDA,
-    ggplot2,
     doParallel,
     foreach,
     parallel,
@@ -26,9 +26,16 @@ Suggests:
     knitr,
     testthat,
     rgl,
-    rmarkdown
+    ggplot2,
+    rmarkdown,
+    SparseM,
+    pliman
 VignetteBuilder:
 		  knitr
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
+Depends: 
+    R (>= 2.10)
+LazyData: true
+LazyDataCompression: xz

R/m_list.R

@@ -0,0 +1,15 @@
+#' Microcebus Teeth as objects for ashapesampler
+#' 
+#' This data set is of 8 of Microcebus mandibular molars characterized as vertices
+#' and the resulting simplicial complex, as used in
+#' Winn-Nunez, et. al. (2023). These were imported from off files with 
+#' function readOFF from the ashapesampler package.
+#'
+#' @format
+#' A named list containing Microcebus molars as mesh3d objects of the Rvcg package 
+#' \describe{
+#'   \item{Vertices}{3 x n matrix containing n vertices as Euclidean coordinates}
+#'   \item{cmplx}{list of vertices, edges, and faces forming simplicial complex}
+#' }
+#' @source <https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/K9A0EG&faces-redirect=true>
+"m_list"

R/m_mesh.R

@@ -0,0 +1,16 @@
+#' Microcebus Teeth as mesh3d objects
+#'
+#' This data set is of 8 mesh3d representations of Microcebus mandibular molars
+#' use in Winn-Nunez, et. al. (2023). These were imported from off files.
+#' See <https://CRAN.R-project.org/package=Rvcg>
+#' for more information about the mesh3d object and the vcgImport function.
+#'
+#' @format
+#' A named list containing Microcebus molars as mesh3d objects of the Rvcg package
+#' \describe{
+#'   \item{vb}{4 x n matrix containing n vertices as homolougous coordinates}
+#'   \item{it}{3 x m matrix containing vertex indices forming triangular faces}
+#'   \item{normals}{4 x n matrix containing vertex normals (homologous coordinates)}
+#' }
+#' @source <https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/K9A0EG&faces-redirect=true>
+"m_mesh"

R/mcmc.R

@@ -1,4 +1,5 @@
 #mcmc sampling
+globalVariables(c("i"))
 
 #' Generate 3D alpha shape
 #'
@@ -20,80 +21,97 @@
 #' @importFrom stats runif
 #' @import doParallel
 #' @import foreach
-generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
-                              afixed = TRUE, mu=NULL, sig = NULL, k_min=3, eps=1e-4,
-                              cores=1){
+generate_ashape3d <- function(point_cloud,
+                              J,
+                              tau,
+                              delta = 0.05,
+                              afixed = TRUE,
+                              mu = NULL,
+                              sig = NULL,
+                              k_min = 3,
+                              eps = 1e-4,
+                              cores = 1) {
   ### Determine the number of Cores for Parallelization ###
-   if(cores > 1){
-     if(cores>parallel::detectCores()){
-       warning("The number of cores you're setting is larger than available cores!")
-       cores <- max(1L, parallel::detectCores(), na.rm = TRUE)}
-   }
-  registerDoParallel(cores=cores)
+  if (cores > 1) {
+    if (cores > parallel::detectCores()) {
+      warning("The number of cores you're setting is larger than available cores!")
+      cores <- max(1L, parallel::detectCores(), na.rm = TRUE)
+    }
+  }
+  registerDoParallel(cores = cores)
   #Check: 3 columns on vertex list
-  if(dim(point_cloud)[2]!=3){
+  if (dim(point_cloud)[2] != 3) {
     stop("Point cloud does not have correct number of columns.")
   }
   n_vert = dim(point_cloud)[1]
-  if(J<=0 || floor(J) !=J){
+  if (J <= 0 || floor(J) != J) {
     stop("J must be positive integer.")
   }
-  if(tau<=0){
+  if (tau <= 0) {
     stop("Tau must be positive real number.")
   }
   #Sample alpha
   my_alpha <- 0
-  if(afixed==FALSE){
-    if(is.null(mu)){
-      mu=tau/3
+  if (afixed == FALSE) {
+    if (is.null(mu)) {
+      mu = tau / 3
     }
-    if(is.null(sig)){
-      sig=tau/12
-    } else if (sig <0){
+    if (is.null(sig)) {
+      sig = tau / 12
+    } else if (sig < 0) {
       stop("sig must be nonnegative value.")
     }
-    if(mu>tau/2 || mu<0){
+    if (mu > tau / 2 || mu < 0) {
       warning("Mean of alpha outside of truncated distribution range for alpha")
     }
-    my_alpha <- truncnorm::rtruncnorm(1, a=0, b=tau/2, mean=mu, sd=sig)
+
+    my_alpha <-
+      truncnorm::rtruncnorm(
+        1,
+        a = 0,
+        b = tau / 2,
+        mean = mu,
+        sd = sig
+      )
   } else {
-    my_alpha <- tau/2-eps
+    my_alpha <- tau / 2 - eps
   }
   #Sample and reject points
-  my_points = matrix(NA, nrow=0, ncol=3)
-  m = n_bound_homology_3D((4/3)*pi*(tau/8)^3, epsilon = my_alpha, tau=tau)
+  my_points = matrix(NA, nrow = 0, ncol = 3)
+  m = n_bound_homology_3D((4 / 3) * pi * (tau / 8) ^ 3, epsilon = my_alpha, tau =
+                            tau)
   my_points = foreach(
-    i = 1:dim(point_cloud)[1],
+    i = 1:n_vert,
     .combine = rbind,
     .export = c("runif_ball_3D", "euclid_dists_point_cloud_3D")
   ) %dopar% {
-  #for (i in 1:n_vert){
-    new_points = runif_ball_3D(m, tau/8)+ cbind(rep(point_cloud[i,1],m),
-                                                rep(point_cloud[i,2],m),
-                                                rep(point_cloud[i,3],m))
-    keep_pts = matrix(NA, nrow=0, ncol=3)
-    for (j in 1:m){
-      dist_list = euclid_dists_point_cloud_3D(new_points[j,], point_cloud)
-      dist_near = dist_list[dist_list < tau/4]
+    #for (i in 1:n_vert){
+    new_points = runif_ball_3D(m, tau / 8) + cbind(rep(point_cloud[i, 1], m),
+                                                   rep(point_cloud[i, 2], m),
+                                                   rep(point_cloud[i, 3], m))
+    keep_pts = matrix(NA, nrow = 0, ncol = 3)
+    for (j in 1:m) {
+      dist_list = euclid_dists_point_cloud_3D(new_points[j, ], point_cloud)
+      dist_near = dist_list[dist_list < tau / 4]
       knn = length(dist_near)
-      if (knn >= k_min*J){
-        keep_pts = rbind(keep_pts, new_points[j,])
+      if (knn >= k_min * J) {
+        keep_pts = rbind(keep_pts, new_points[j, ])
       } else if (knn >= k_min) {
-        a_prob = 1-exp(-(knn-k_min)*2/(J*k_min))
-        if (runif(1)<a_prob){
-          keep_pts = rbind(keep_pts, new_points[j,])
+        a_prob = 1 - exp(-(knn - k_min) * 2 / (J * k_min))
+        if (runif(1) < a_prob) {
+          keep_pts = rbind(keep_pts, new_points[j, ])
+        }
       }
     }
-    }
     keep_pts
   }
   my_points = unique(my_points) #keeps error free if necessary.
-  if(dim(my_points)[1]<5){
+  if (dim(my_points)[1] < 5) {
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 5.")
   }
-  rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
+  rr = dim(my_points)[1] / (m * dim(point_cloud)[1])
   print(paste0("Acceptance Rate is ", rr))
-  new_ashape <- alphashape3d::ashape3d(my_points, alpha=tau-eps)
+  new_ashape <- alphashape3d::ashape3d(my_points, alpha = tau - eps)
   return(new_ashape)
 }
 
@@ -117,77 +135,95 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
 #' @importFrom stats runif
 #' @import doParallel
 #' @import foreach
-generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
-                              afixed=TRUE, mu=NULL, sig = NULL, k_min=2, eps=1e-4,
-                              cores=1){
+generate_ashape2d <- function(point_cloud,
+                              J,
+                              tau,
+                              delta = 0.05,
+                              afixed = TRUE,
+                              mu = NULL,
+                              sig = NULL,
+                              k_min = 2,
+                              eps = 1e-4,
+                              cores = 1) {
   ### Determine the number of Cores for Parallelization ###
-  if(cores > 1){
-    if(cores>parallel::detectCores()){
+  if (cores > 1) {
+    if (cores > parallel::detectCores()) {
       warning("The number of cores you're setting is larger than available cores!")
-      cores <- max(1L, parallel::detectCores(), na.rm = TRUE)}
+      cores <- max(1L, parallel::detectCores(), na.rm = TRUE)
+    }
   }
-  registerDoParallel(cores=cores)
+  registerDoParallel(cores = cores)
   #Check: 2 columns on vertex list
-  if(dim(point_cloud)[2]!=2){
+  if (dim(point_cloud)[2] != 2) {
     stop("Point cloud does not have correct number of columns.")
   }
   n_vert = dim(point_cloud)[1]
-  if(J<=0 || floor(J) !=J){
+  if (J <= 0 || floor(J) != J) {
     stop("J must be positive integer.")
   }
-  if(tau<=0){
+  if (tau <= 0) {
     stop("Tau must be positive real number.")
   }
   #Sample alpha
-  my_alpha=0
-  if(afixed==FALSE){
-    if(is.null(mu)){
-      mu=tau/3
+  my_alpha = 0
+  if (afixed == FALSE) {
+    if (is.null(mu)) {
+      mu = tau / 3
     }
-    if(is.null(sig)){
-      sig=tau/12
-    } else if (sig <0){
+    if (is.null(sig)) {
+      sig = tau / 12
+    } else if (sig < 0) {
       stop("sig must be nonnegative value.")
     }
-    if(mu>tau/2 || mu<0){
+    if (mu > tau / 2 || mu < 0) {
       warning("Mean of alpha outside of truncated distribution range for alpha")
     }
-    my_alpha <- truncnorm::rtruncnorm(1, a=0, b=tau/2, mean=mu, sd=sig)
+    my_alpha <-
+      truncnorm::rtruncnorm(
+        1,
+        a = 0,
+        b = tau / 2,
+        mean = mu,
+        sd = sig
+      )
   } else {
-    my_alpha <- tau/2-eps
+    my_alpha <- tau / 2 - eps
   }
   #Sample and reject points
-  my_points = matrix(NA, nrow=0, ncol=2)
+  my_points = matrix(NA, nrow = 0, ncol = 2)
   #Initialize by taking point from point cloud.
-  m = n_bound_homology_2D(pi*(tau/8)^2, epsilon = my_alpha, tau=tau)
-
- my_points = foreach(i=1:n_vert, .combine=rbind,
-          .export = c("runif_disk", "euclid_dists_point_cloud_2D"))%dopar%{
+  m = n_bound_homology_2D(pi * (tau / 8) ^ 2, epsilon = my_alpha, tau =
+                            tau)
 
-  #for(i in 1:n_vert){
-    new_points = runif_disk(m, tau/8)+ cbind(rep(point_cloud[i,1],m), rep(point_cloud[i,2],m))
-    keep_pts = matrix(NA, nrow=0, ncol=2)
-     for (j in 1:m){
-       dist_list = euclid_dists_point_cloud_2D(new_points[j,], point_cloud)
-       dist_near = dist_list[dist_list < tau/4]
-       knn = length(dist_near)
-       if (knn >= k_min*J){
-        keep_pts = rbind(keep_pts, new_points[j,])
+  my_points = foreach(
+    i = 1:n_vert,
+    .combine = rbind,
+    .export = c("runif_disk", "euclid_dists_point_cloud_2D")
+  ) %dopar% {
+    #for(i in 1:n_vert){
+    new_points = runif_disk(m, tau / 8) + cbind(rep(point_cloud[i, 1], m), rep(point_cloud[i, 2], m))
+    keep_pts = matrix(NA, nrow = 0, ncol = 2)
+    for (j in 1:m) {
+      dist_list = euclid_dists_point_cloud_2D(new_points[j, ], point_cloud)
+      dist_near = dist_list[dist_list < tau / 4]
+      knn = length(dist_near)
+      if (knn >= k_min * J) {
+        keep_pts = rbind(keep_pts, new_points[j, ])
       } else if (knn >= k_min) {
-        a_prob = 1-exp(-(knn-k_min)*2/(J*k_min))
-        if (runif(1)<a_prob){
-          keep_pts = rbind(keep_pts, new_points[j,])
+        a_prob = 1 - exp(-(knn - k_min) * 2 / (J * k_min))
+        if (runif(1) < a_prob) {
+          keep_pts = rbind(keep_pts, new_points[j, ])
         }
-     }
-     }
+      }
+    }
     keep_pts
-          }
+  }
   my_points = unique(my_points) #keeps error free if necessary.
-  if(dim(my_points)[1]<3){
+  if (dim(my_points)[1] < 3) {
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 3.")
   }
-  rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
+  rr = dim(my_points)[1] / (m * dim(point_cloud)[1])
   print(paste0("Acceptance Rate is ", rr))
-  new_ashape <- alphahull::ashape(my_points, alpha=tau-eps)
+  new_ashape <- alphahull::ashape(my_points, alpha = tau - eps)
   return(new_ashape)
 }

R/t_list.R

@@ -0,0 +1,15 @@
+#' Tarsius Teeth as objects for ashapesampler
+#' 
+#' This data set is of 7 of Tarsius mandibular molars characterized as vertices
+#' and the resulting simplicial complex, as used in
+#' Winn-Nunez, et. al. (2023). These were imported from off files with 
+#' function readOFF from the ashapesampler package.
+#'
+#' @format
+#' A named list containing Microcebus molars as mesh3d objects of the Rvcg package 
+#' \describe{
+#'   \item{Vertices}{3 x n matrix containing n vertices as Euclidean coordinates}
+#'   \item{cmplx}{list of vertices, edges, and faces forming simplicial complex}
+#' }
+#' @source <https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/K9A0EG&faces-redirect=true>
+"t_list"