Merge pull request #22 from lcrawlab/new_functions

editing mcmc functions

R/mcmc.R

@@ -12,6 +12,10 @@ globalVariables(c("i"))
 #' @param mu mean of truncated distribution from which alpha sampled; default tau/3
 #' @param sig standard deviation of truncated distribution from which alpha
 #'              sampled; default tau/12
+#' @param sample_rad radius of ball around each point in point cloud from which to
+#'                   sample; default tau/8
+#' @param acc_rad radius of ball to check around potential sampled points for whether
+#'                to accept or reject new point; default tau/4
 #' @param k_min number of points needed in radius 2 alpha of point cloud to accept a sample
 #' @param eps amount to subtract from tau/2 to give alpha. Defaul 1e-4.
 #' @param cores number of cores for parallelizing. Default 1.
@@ -28,6 +32,8 @@ generate_ashape3d <- function(point_cloud,
                               afixed = TRUE,
                               mu = NULL,
                               sig = NULL,
+                              sample_rad=NULL,
+                              acc_rad = NULL,
                               k_min = 3,
                               eps = 1e-4,
                               cores = 1) {
@@ -76,23 +82,35 @@ generate_ashape3d <- function(point_cloud,
   } else {
     my_alpha <- tau / 2 - eps
   }
+  #Sampling radius
+  if(is.null(sample_rad)){
+    sample_rad = tau/8
+  } else if (sample_rad <= 0){
+    stop("Sampling radius must be larger than 0.")
+  }
+  #acceptance radius
+  if(is.null(acc_rad)){
+    acc_rad = tau/4
+  } else if (acc_rad <= 0){
+    stop("Acceptance radius must be larger than 0.")
+  }
   #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 =
+  m = n_bound_homology_3D((4 / 3) * pi * (sample_rad) ^ 3, epsilon = my_alpha, tau =
                             tau)
   my_points = foreach(
     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),
+    new_points = runif_ball_3D(m, sample_rad) + 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]
+      dist_near = dist_list[dist_list < acc_rad]
       knn = length(dist_near)
       if (knn >= k_min * J) {
         keep_pts = rbind(keep_pts, new_points[j, ])
@@ -107,7 +125,7 @@ generate_ashape3d <- function(point_cloud,
   }
   my_points = unique(my_points) #keeps error free if necessary.
   if (dim(my_points)[1] < 5) {
-    stop("Not enough points accepted to make a shape. Need at least 5. Check tau and k_min parameters to 
+    stop("Not enough points accepted to make a shape. Need at least 5. Check tau and k_min parameters to
          increase probability of acceptance.")
   }
   rr = dim(my_points)[1] / (m * dim(point_cloud)[1])
@@ -127,6 +145,10 @@ generate_ashape3d <- function(point_cloud,
 #' @param mu mean of truncated distribution from which alpha sampled; default tau/3
 #' @param sig standard deviation of truncated distribution from which alpha
 #'              sampled; default tau/12
+#' @param sample_rad radius of ball around each point in point cloud from which to
+#'                   sample; default tau/8
+#' @param acc_rad radius of ball to check around potential sampled points for whether
+#'                to accept or reject new point; default tau/4
 #' @param k_min number of points needed in radius tau of point cloud to accept a sample
 #' @param eps amount to subtract from tau/2 to give alpha. Defaul 1e-4.
 #' @param cores number of computer cores for parallelizing. Default 1.
@@ -143,6 +165,8 @@ generate_ashape2d <- function(point_cloud,
                               afixed = TRUE,
                               mu = NULL,
                               sig = NULL,
+                              sample_rad=NULL,
+                              acc_rad = NULL,
                               k_min = 2,
                               eps = 1e-4,
                               cores = 1) {
@@ -190,10 +214,22 @@ generate_ashape2d <- function(point_cloud,
   } else {
     my_alpha <- tau / 2 - eps
   }
+  #Sampling radius
+  if(is.null(sample_rad)){
+    sample_rad = tau/8
+  } else if (sample_rad <=0){
+    stop("Sampling radius must be larger than 0.")
+  }
+  #acceptance radius
+  if(is.null(acc_rad)){
+    acc_rad = tau/4
+  } else if (acc_rad <= 0){
+    stop("Acceptance radius must be larger than 0.")
+  }
   #Sample and reject points
   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 =
+  m = n_bound_homology_2D(pi * (sample_rad) ^ 2, epsilon = my_alpha, tau =
                             tau)
 
   my_points = foreach(
@@ -202,11 +238,11 @@ generate_ashape2d <- function(point_cloud,
     .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))
+    new_points = runif_disk(m, sample_rad) + 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]
+      dist_near = dist_list[dist_list < acc_rad]
       knn = length(dist_near)
       if (knn >= k_min * J) {
         keep_pts = rbind(keep_pts, new_points[j, ])
@@ -221,7 +257,7 @@ generate_ashape2d <- function(point_cloud,
   }
   my_points = unique(my_points) #keeps error free if necessary.
   if (dim(my_points)[1] < 3) {
-    stop("Not enough points accepted to make a shape. Need at least 3. Check tau and k_min parameters to 
+    stop("Not enough points accepted to make a shape. Need at least 3. Check tau and k_min parameters to
          increase probability of acceptance.")
   }
   rr = dim(my_points)[1] / (m * dim(point_cloud)[1])

tests/testthat/test_mcmc.R

@@ -21,6 +21,10 @@ test_that("mcmc errors", {
   #expect_warning(generate_ashape3d(points3, J=2, tau=5, cores=120))
   expect_error(generate_ashape2d(points2, J=2, tau=5, k_min=100))
   expect_error(generate_ashape3d(points3, J=2, tau=6, k_min=1000))
+  expect_error(generate_ashape2d(points2, N, tau, sample_rad=-1))
+  expect_error(generate_ashape2d(points2, N, tau, acc_rad=-1))
+  expect_error(generate_ashape3d(points3, N, tau, sample_rad=-1))
+  expect_error(generate_ashape3d(points3, N, tau, acc_rad=-1))
 })
 
 test_that("2D mcmc smooth runs", {