diff --git a/R/2D.R b/R/2D.R
index 8b8ed59..b815527 100644
--- a/R/2D.R
+++ b/R/2D.R
@@ -184,7 +184,7 @@ euclid_dists_point_cloud_2D <- function(point, point_cloud){
   }
   dist_vec = vector("numeric",m)
   for (j in 1:m){
-    sqr_dist = (point[1]-point_cloud[j,1])^2+(point[2]-point_cloud[j,2])^2
+    sqr_dist = as.numeric((point[1]-point_cloud[j,1])^2+(point[2]-point_cloud[j,2])^2)
     dist_vec[j] = sqrt(sqr_dist)
   }
   return(dist_vec)
diff --git a/R/mcmc.R b/R/mcmc.R
index 61036a1..eeca19e 100644
--- a/R/mcmc.R
+++ b/R/mcmc.R
@@ -68,7 +68,9 @@ generate_ashape3d <- function(point_cloud, J, tau, delta=0.05,
     .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)+rep(point_cloud[i,], each=m)
+    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)
@@ -163,7 +165,7 @@ generate_ashape2d <- function(point_cloud, J, tau, delta=0.05,
           .export = c("runif_disk", "euclid_dists_point_cloud_2D"))%dopar%{
 
   #for(i in 1:n_vert){
-    new_points = runif_disk(m, tau/8)+rep(point_cloud[i,], each=m)
+    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)
diff --git a/R/tau_bound.R b/R/tau_bound.R
index db9bf82..02833bc 100644
--- a/R/tau_bound.R
+++ b/R/tau_bound.R
@@ -68,16 +68,16 @@ tau_bound <- function(v_list, complex, extremes=NULL, cores = 1){
     m=length(extremes)
   }
   dist_matrix = as.matrix(stats::dist(v_list))
-  e_list = extract_complex_edges(complex,m)
+  e_list = extract_complex_edges(complex,n)
   if(is.null(e_list)){
     return(min(dist_matrix[dist_matrix>0]))
   }
-  f_list = extract_complex_faces(complex,m)
+  f_list = extract_complex_faces(complex,n)
   f_circ = circumcenter_face(v_list, f_list)
   t_list = NULL
   t_circ = NULL
   if(dimension>2){
-    t_list = extract_complex_tet(complex,m)
+    t_list = extract_complex_tet(complex,n)
     t_circ = circumcenter_tet(v_list, t_list)
   }
   tau_vec=vector("numeric", m)
@@ -108,16 +108,16 @@ tau_bound <- function(v_list, complex, extremes=NULL, cores = 1){
       dist_vec = dist_vec_point[edge_list_zoom]
       dist_vec_b = c()
       if (dimension == 2){
-        if(!is.null(face_list_zoom)){
+        if(length(face_list_zoom)>0){
           points = matrix(f_circ[face_list_zoom,], ncol=2)
           dist_vec_b = c(dist_vec_b, 2*euclid_dists_point_cloud_2D(v_list[i,],
                                                               points ))
         }
       } else {
-        if(!is.null(face_list_zoom)){
+        if(length(face_list_zoom)>0){
           points = matrix(f_circ[face_list_zoom,], ncol=3)
           dist_vec_b = 2*euclid_dists_point_cloud_3D(v_list[i,],points)
-          if(!is.null(tet_list_zoom)){
+          if(length(tet_list_zoom)>0){
           points = matrix(t_circ[tet_list_zoom,], ncol=3)
           dist_vec_b = c(dist_vec_b, 2*euclid_dists_point_cloud_3D(v_list[i,], points))
           }