Merge pull request #17 from lcrawlab/new_functions

finalizing algorithm

DESCRIPTION

@@ -20,11 +20,11 @@ Imports:
     ggplot2,
     doParallel,
     foreach,
-    parallel
+    parallel,
+    dplyr
 Suggests:
     knitr,
     testthat,
-    dplyr,
     rgl,
     rmarkdown
 VignetteBuilder:

NAMESPACE

@@ -8,6 +8,7 @@ export(euclid_dists_point_cloud_3D)
 export(extract_complex_edges)
 export(extract_complex_faces)
 export(extract_complex_tet)
+export(extreme_pts)
 export(generate_ashape2d)
 export(generate_ashape3d)
 export(get_alpha_complex)
@@ -31,5 +32,6 @@ import(doParallel)
 import(foreach)
 importFrom(Rvcg,vcgGetEdge)
 importFrom(Rvcg,vcgImport)
+importFrom(dplyr,setdiff)
 importFrom(stats,na.omit)
 importFrom(stats,runif)

R/2D.R

@@ -23,7 +23,7 @@ circle_overlap_cc <- function(alpha, r = 1) {
 
 #' Circle Overlap Inner Annulus
 #'
-#' Circle overlap ia (inner annulus) calculates area needed to subtract 
+#' Circle overlap ia (inner annulus) calculates area needed to subtract
 #' when calculating area of overlap of annulus and circle.
 #'
 #' @param alpha radius of circle
@@ -39,7 +39,7 @@ circle_overlap_ia <- function(alpha, R, r){
   theta_2 = 2*acos(((R^2+r^2-alpha^2))/(2*r*R))
   a_1 = alpha^2*0.5*(theta_1 - sin(theta_1))
   a_2 = 0
-  if (alpha^2 < R^2 + r^2){ 
+  if (alpha^2 < R^2 + r^2){
     a_2 = r^2*0.5*(theta_2 - sin(theta_2))
   } else {
     a_2 = r^2*0.5*(theta_2 + sin((2*pi-theta_2)))
@@ -51,7 +51,7 @@ circle_overlap_ia <- function(alpha, R, r){
 #'
 #' This function calculates the minimum coverage percentage of an alpha ball over the bounded
 #' area being considered. 0 is no coverage, 1 means complete coverage.
-#' For the square, r is the length of the side. For circle, r is the radius. For 
+#' For the square, r is the length of the side. For circle, r is the radius. For
 #' the annulus, r and min_r are the two radii.
 #' @param alpha radius of alpha ball
 #' @param r length of square, radius of circle, or outer radius of annulus
@@ -101,7 +101,7 @@ calc_overlap_2D <- function(alpha, r=1, rmin=0.01, bound="square"){
       return(1)
     }
   } else {
-    stop("Not a valid bound for the manifold. Please enter 'square', 'circle', or 
+    stop("Not a valid bound for the manifold. Please enter 'square', 'circle', or
           'annulus'. Default is 'square'.")
   }
 }
@@ -133,13 +133,13 @@ n_bound_connect_2D <- function(alpha, delta=0.05, r=1, rmin=0.01, bound="square"
 
 #from Niyogi et al 2008
 #' n Bound Homology 2D
-#' 
-#' #' Function returns the minimum number of points to preserve the homology with 
+#'
+#' #' Function returns the minimum number of points to preserve the homology with
 #' an open cover of radius alpha.
 #'
 #' @param area area of manifold from which points being sampled
 #' @param epsilon size of balls of cover
-#' @param tau number bound 
+#' @param tau number bound
 #' @param delta probability of not recovering homology
 #'
 #' @return n, number of points needed
@@ -222,20 +222,20 @@ runif_square <- function(n, xmin=0, xmax=1, ymin=0, ymax=1){
 }
 
 #' Uniform sampling from disk
-#' 
+#'
 #' Returns points uniformly sampled from disk of radius r in plane
 #'
 #' @param n number of points to sample
 #' @param r radius of disk
 #'
 #' @return points n by 2 matrix of points sampled
 #'
-#' @examples 
+#' @examples
 #' # Sample 100 points from unit disk
 #' runif_disk(100)
 #' # Sample 100 points from disk of radius 0.7
 #' runif_disk(100, 0.7)
-#' @export 
+#' @export
 runif_disk <- function(n, r=1){
   if(n<=0 || floor(n) !=n || r<=0){
     stop("n must be positive integer, and r must be a positive real number.")
@@ -249,7 +249,7 @@ runif_disk <- function(n, r=1){
 }
 
 #' Uniform Sampling from Annulus
-#' 
+#'
 #' Returns points uniformly sampled from annulus in plane
 #'
 #' @param n number of points to sample

R/complex.R

@@ -47,7 +47,8 @@ extract_complex_edges <- function(complex, n_vert=0){
   #First n_vert entries are the vertices, no need to check
   for (k in (n_vert+1):m){
     if(length(as.vector(complex[[k]]))==2){
-      edge_list[(k-n_vert),] <- complex[[k]]
+      temp <- sort(complex[[k]])
+      edge_list[(k-n_vert),] <- temp
     }
   }
   edge_list = as.data.frame(na.omit(edge_list))
@@ -82,7 +83,8 @@ extract_complex_faces <- function(complex, n_vert=0){
   #First n_vert entries are the vertices, no need to check
   for (k in (n_vert+1):m){
     if(length(as.vector(complex[[k]]))==3){
-      face_list[(k-n_vert),] <- complex[[k]]
+      temp <- sort(complex[[k]])
+      face_list[(k-n_vert),] <- temp
     }
   }
   face_list = as.data.frame(na.omit(face_list))
@@ -115,7 +117,8 @@ extract_complex_tet <- function(complex, n_vert=0){
   #First n_vert entries are the vertices, no need to check
   for (k in (n_vert+1):m){
     if(length(as.vector(complex[[k]]))==4){
-      tet_list[(k-n_vert),] <- complex[[k]]
+      temp <- sort(complex[[k]])
+      tet_list[(k-n_vert),] <- temp
     }
   }
   tet_list = as.data.frame(na.omit(tet_list))

R/mcmc.R

@@ -5,7 +5,7 @@
 #' @param point_cloud 3 column matrix of all points from all shapes in initial
 #'                    data set
 #' @param J number of shapes in initial data set
-#' @param tau tau bound
+#' @param tau tau bound for the shapes
 #' @param delta probability of not preserving homology; default is 0.05
 #' @param afixed boolean, whether to sample alpha or leave fixed based on tau. Default FALSE
 #' @param mu mean of truncated distribution from which alpha sampled; default tau/3
@@ -59,7 +59,6 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
   } else {
     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)
@@ -73,7 +72,7 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
     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]
+      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,])
@@ -90,7 +89,9 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
   if(dim(my_points)[1]<5){
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 5.")
   }
-  new_ashape <- alphashape3d::ashape3d(my_points, alpha=tau-eps)
+  rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
+  print(paste0("Acceptance Rate is ", rr))
+  new_ashape <- alphashape3d::ashape3d(my_points, alpha=my_alpha)
   return(new_ashape)
 }
 
@@ -99,7 +100,7 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
 #' @param point_cloud 2 column matrix of all points from all shapes in initial
 #'                    data set
 #' @param J number of shapes in initial (sub) data set
-#' @param tau tau bound
+#' @param tau tau bound vector for shapes input
 #' @param delta probability of not preserving homology; default is 0.05
 #' @param afixed boolean, whether to sample alpha or leave fixed based on tau. Default FALSE
 #' @param mu mean of truncated distribution from which alpha sampled; default tau/3
@@ -166,7 +167,7 @@ generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
     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]
+       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,])
@@ -182,6 +183,8 @@ generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
   if(dim(my_points)[1]<3){
     stop("Not enough points accepted in MCMC walk to make a shape. Need at least 3.")
   }
-  new_ashape <- alphahull::ashape(my_points, alpha=tau-eps)
+  rr = dim(my_points)[1]/(m*dim(point_cloud)[1])
+  print(paste0("Acceptance Rate is ", rr))
+  new_ashape <- alphahull::ashape(my_points, alpha=my_alpha)
   return(new_ashape)
 }