diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION
new file mode 100644
index 0000000..958d534
--- /dev/null
+++ b/CRAN-SUBMISSION
@@ -0,0 +1,3 @@
+Version: 0.4.1
+Date: 2022-11-22 06:42:52 UTC
+SHA: 9196928bc4eda31af1e25c30dc6e71912b87cf58
diff --git a/DESCRIPTION b/DESCRIPTION
index 843e41c..763e8c1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: edgebundle
 Title: Algorithms for Bundling Edges in Networks and Visualizing Flow and Metro Maps
-Version: 0.4.1
+Version: 0.4.2
 Authors@R: 
     person(given = "David",
            family = "Schoch",
@@ -19,7 +19,7 @@ Config/testthat/edition: 2
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 LinkingTo: 
     Rcpp
 Imports: 
diff --git a/NEWS.md b/NEWS.md
index a67deee..f9090cb 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# edgebundle 0.4.2
+
+* clean up codebase
+
 # edgebundle 0.4.1
 
 * added package level docs
diff --git a/R/bundle_force.R b/R/bundle_force.R
index 8268546..9ea25ec 100644
--- a/R/bundle_force.R
+++ b/R/bundle_force.R
@@ -20,39 +20,49 @@
 #' Holten, Danny, and Jarke J. Van Wijk. "Force-Directed Edge Bundling for Graph Visualization." Computer Graphics Forum (Blackwell Publishing Ltd) 28, no. 3 (2009): 983-990.
 #' @examples
 #' library(igraph)
-#' g  <- graph_from_edgelist(matrix(c(1,12,2,11,3,10,4,9,5,8,6,7),ncol = 2,byrow = TRUE),FALSE)
-#' xy <- cbind(c(rep(0,6),rep(1,6)),c(1:6,1:6))
-#' edge_bundle_force(g,xy)
+#' g <- graph_from_edgelist(
+#'     matrix(c(
+#'         1, 12, 2, 11, 3, 10,
+#'         4, 9, 5, 8, 6, 7
+#'     ), ncol = 2, byrow = TRUE), FALSE
+#' )
+#' xy <- cbind(c(rep(0, 6), rep(1, 6)), c(1:6, 1:6))
+#' edge_bundle_force(g, xy)
 #' @export
 
-edge_bundle_force <- function(object,xy, K=1,C=6,P=1,S=0.04,
-                              P_rate=2,I=50,I_rate=2/3,
+edge_bundle_force <- function(object, xy, K = 1, C = 6, P = 1, S = 0.04,
+                              P_rate = 2, I = 50, I_rate = 2 / 3,
                               compatibility_threshold = 0.6,
                               eps = 1e-8) {
-  #initialize matrix with coordinates
-  edges_xy <- convert_edges(object,xy)
-  m <- nrow(edges_xy)
+    # initialize matrix with coordinates
+    edges_xy <- convert_edges(object, xy)
+    m <- nrow(edges_xy)
 
-  #initialize edge subdivision list
-  elist <- unname(lapply(split(edges_xy, rep(1:nrow(edges_xy), ncol(edges_xy))),
-                         function(y) matrix(y,2,2,byrow = TRUE)))
+    # initialize edge subdivision list
+    elist <- unname(lapply(
+        split(edges_xy, rep(seq_len(nrow(edges_xy)), ncol(edges_xy))),
+        function(y) matrix(y, 2, 2, byrow = TRUE)
+    ))
 
-  #main force bundling routine
-  elist <- force_bundle_iter(edges_xy,elist,K,C,P,P_rate,
-                             S,I, I_rate,compatibility_threshold, eps)
+    # main force bundling routine
+    elist <- force_bundle_iter(
+        edges_xy, elist, K, C, P, P_rate,
+        S, I, I_rate, compatibility_threshold, eps
+    )
 
-  # assemble data frame
-  segments <- nrow(elist[[1]])
+    # assemble data frame
+    segments <- nrow(elist[[1]])
 
-  idx <- seq(0, 1, length.out = segments)
-  data_bundle <- as.data.frame(cbind(
-    do.call("rbind",elist),
-    rep(idx,m),
-    rep(1:m,each=segments)))
+    idx <- seq(0, 1, length.out = segments)
+    data_bundle <- as.data.frame(cbind(
+        do.call("rbind", elist),
+        rep(idx, m),
+        rep(1:m, each = segments)
+    ))
 
-  names(data_bundle) <- c("x","y","index","group")
+    names(data_bundle) <- c("x", "y", "index", "group")
 
-  data_bundle
+    data_bundle
 }
 
 #' @importFrom Rcpp sourceCpp
diff --git a/R/bundle_hammer.R b/R/bundle_hammer.R
index 16d5054..69794d2 100644
--- a/R/bundle_hammer.R
+++ b/R/bundle_hammer.R
@@ -14,40 +14,38 @@
 #' @seealso [edge_bundle_force],[edge_bundle_stub], [edge_bundle_path]
 #' @export
 
-edge_bundle_hammer <- function(object,xy,bw=0.05,decay=0.7){
-  if (!requireNamespace('reticulate', quietly = TRUE)) {
-    stop('The `reticulate` package is required for this functionality')
-  }
-  if(any(class(object)=="igraph")){
-    if (!requireNamespace('igraph', quietly = TRUE)) {
-      stop('The `igraph` package is required for this functionality')
+edge_bundle_hammer <- function(object, xy, bw = 0.05, decay = 0.7) {
+    if (!requireNamespace("reticulate", quietly = TRUE)) {
+        stop("The `reticulate` package is required for this functionality")
     }
-    nodes <- data.frame(name=paste0("node",0:(igraph::vcount(object)-1)),x=xy[,1],y=xy[,2])
-    el <- igraph::get.edgelist(object,names = FALSE)
-    el1 <- data.frame(source=el[,1]-1,target=el[,2]-1)
-
-  } else if(any(class(object)=="tbl_graph")){
-    if (!requireNamespace('tidygraph', quietly = TRUE)) {
-      stop('The `tidygraph` package is required for this functionality')
+    if (any(class(object) == "igraph")) {
+        if (!requireNamespace("igraph", quietly = TRUE)) {
+            stop("The `igraph` package is required for this functionality")
+        }
+        nodes <- data.frame(name = paste0("node", 0:(igraph::vcount(object) - 1)), x = xy[, 1], y = xy[, 2])
+        el <- igraph::get.edgelist(object, names = FALSE)
+        el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1)
+    } else if (any(class(object) == "tbl_graph")) {
+        if (!requireNamespace("tidygraph", quietly = TRUE)) {
+            stop("The `tidygraph` package is required for this functionality")
+        }
+        object <- tidygraph::as.igraph(object)
+        nodes <- data.frame(name = paste0("node", 0:(igraph::vcount(object) - 1)), x = xy[, 1], y = xy[, 2])
+        el <- igraph::get.edgelist(object, names = FALSE)
+        el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1)
+    } else if (any(class(object) == "network")) {
+        nodes <- data.frame(name = paste0("node", 0:(network::get.network.attribute(object, "n") - 1)), x = xy[, 1], y = xy[, 2])
+        el <- network::as.edgelist(object)
+        el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1)
+    } else {
+        stop("only `igraph`, `network` or `tbl_graph` objects supported.")
     }
-    object <- tidygraph::as.igraph(object)
-    nodes <- data.frame(name=paste0("node",0:(igraph::vcount(object)-1)),x=xy[,1],y=xy[,2])
-    el <- igraph::get.edgelist(object,names = FALSE)
-    el1 <- data.frame(source=el[,1]-1,target=el[,2]-1)
-
-  } else if(any(class(object)=="network")){
-    nodes <- data.frame(name=paste0("node",0:(network::get.network.attribute(object,"n")-1)),x=xy[,1],y=xy[,2])
-    el <- network::as.edgelist(object)
-    el1 <- data.frame(source=el[,1]-1,target=el[,2]-1)
-  } else{
-    stop("only `igraph`, `network` or `tbl_graph` objects supported.")
-  }
-  data_bundle <- shader_env$datashader_bundling$hammer_bundle(nodes,el1,initial_bandwidth = bw,decay = decay)
-  data_bundle$group <- is.na(data_bundle$y)+0
-  data_bundle$group <- cumsum(data_bundle$group)+1
-  data_bundle <- data_bundle[!is.na(data_bundle$y),]
-  data_bundle$index <- unlist(sapply(table(data_bundle$group),function(x) seq(0,1,length.out=x)))
-  data_bundle[,c("x","y","index","group")]
+    data_bundle <- shader_env$datashader_bundling$hammer_bundle(nodes, el1, initial_bandwidth = bw, decay = decay)
+    data_bundle$group <- is.na(data_bundle$y) + 0
+    data_bundle$group <- cumsum(data_bundle$group) + 1
+    data_bundle <- data_bundle[!is.na(data_bundle$y), ]
+    data_bundle$index <- unlist(sapply(table(data_bundle$group), function(x) seq(0, 1, length.out = x)))
+    data_bundle[, c("x", "y", "index", "group")]
 }
 
 #' @title install python dependencies for hammer bundling
@@ -60,17 +58,17 @@ edge_bundle_hammer <- function(object,xy,bw=0.05,decay=0.7){
 #' @export
 #'
 install_bundle_py <- function(method = "auto", conda = "auto") {
-  if (!requireNamespace('reticulate', quietly = TRUE)) {
-    stop('The `reticulate` package is required for this functionality')
-  }
-  reticulate::py_install("datashader", method = method, conda = conda, pip = TRUE)
-  reticulate::py_install("scikit-image", method = method, conda = conda, pip = TRUE)
+    if (!requireNamespace("reticulate", quietly = TRUE)) {
+        stop("The `reticulate` package is required for this functionality")
+    }
+    reticulate::py_install("datashader", method = method, conda = conda, pip = TRUE)
+    reticulate::py_install("scikit-image", method = method, conda = conda, pip = TRUE)
 }
 
 # Environment for globals
 shader_env <- new.env(parent = emptyenv())
 
 .onLoad <- function(libname, pkgname) {
-  reticulate::configure_environment(pkgname)
-  assign("datashader_bundling", reticulate::import("datashader.bundling", delay_load = TRUE), shader_env)
+    reticulate::configure_environment(pkgname)
+    assign("datashader_bundling", reticulate::import("datashader.bundling", delay_load = TRUE), shader_env)
 }
diff --git a/R/bundle_path.R b/R/bundle_path.R
index 617cbfa..16cfd84 100644
--- a/R/bundle_path.R
+++ b/R/bundle_path.R
@@ -14,93 +14,99 @@
 #' Wallinger, M., Archambault, D., Auber, D., Nollenburg, M., & Peltonen, J. (2021). Edge-Path Bundling: A Less Ambiguous Edge Bundling Approach. IEEE Transactions on Visualization and Computer Graphics.
 #' @examples
 #' library(igraph)
-#' g  <- graph_from_edgelist(matrix(c(1,2,1,6,1,4,2,3,3,4,4,5,5,6),ncol = 2,byrow = TRUE),FALSE)
-#' xy <- cbind(c(0,10,25,40,50,50),c(0,15,25,15,0,-10))
-#' edge_bundle_path(g,xy)
+#' g <- graph_from_edgelist(matrix(c(
+#'     1, 2, 1, 6,
+#'     1, 4, 2, 3, 3, 4, 4, 5, 5, 6
+#' ), ncol = 2, byrow = TRUE), FALSE)
+#' xy <- cbind(c(0, 10, 25, 40, 50, 50), c(0, 15, 25, 15, 0, -10))
+#' edge_bundle_path(g, xy)
 #' @export
 
-edge_bundle_path <- function(g,xy,max_distortion = 2,weight_fac = 2,segments = 20){
-  # preprocess
-  if(!igraph::is.igraph(g)){
-    stop("edge_bundle_path requires the input graph to be an ingraph object")
-  }
-  m <- igraph::ecount(g)
-  lock <- rep(FALSE,m)
-  skip <- rep(FALSE,m)
-
-  el <- igraph::get.edgelist(g,names = FALSE)
-  exy <- cbind(xy[el[,1],1],xy[el[,1],2],
-               xy[el[,2],1],xy[el[,2],2])
-  elen <- sqrt((exy[,1]-exy[,3])^2+(exy[,2]-exy[,4])^2)
-  weights <- elen^weight_fac
-  sortedEdges <- order(weights,decreasing = TRUE)
-  cpoints <- vector("list",m)
-  #iterate
-  for(e in sortedEdges){
-    s <- el[e,1]
-    t <- el[e,2]
-    cpoints[[e]] <- xy[c(s,t),]
-    if(lock[e]){
-      next()
-    }
-    skip[e] <- TRUE
-    g1 <- igraph::delete.edges(g,which(skip))
-    sp_verts <- suppressWarnings(igraph::shortest_paths(g1,s,t,weights = weights[!skip])$vpath[[1]])
-    if(length(sp_verts)<2){
-      skip[e] <- FALSE
-      next
+edge_bundle_path <- function(g, xy, max_distortion = 2, weight_fac = 2, segments = 20) {
+    # preprocess
+    if (!igraph::is.igraph(g)) {
+        stop("edge_bundle_path requires the input graph to be an ingraph object")
     }
-    sp_len <- path_length(sp_verts,xy)
-    if(sp_len >= max_distortion * elen[e]){
-      skip[e] <- FALSE
-      next
+    m <- igraph::ecount(g)
+    lock <- rep(FALSE, m)
+    skip <- rep(FALSE, m)
+
+    el <- igraph::get.edgelist(g, names = FALSE)
+    exy <- cbind(
+        xy[el[, 1], 1], xy[el[, 1], 2],
+        xy[el[, 2], 1], xy[el[, 2], 2]
+    )
+    elen <- sqrt((exy[, 1] - exy[, 3])^2 + (exy[, 2] - exy[, 4])^2)
+    weights <- elen^weight_fac
+    sortedEdges <- order(weights, decreasing = TRUE)
+    cpoints <- vector("list", m)
+    # iterate
+    for (e in sortedEdges) {
+        s <- el[e, 1]
+        t <- el[e, 2]
+        cpoints[[e]] <- xy[c(s, t), ]
+        if (lock[e]) {
+            next()
+        }
+        skip[e] <- TRUE
+        g1 <- igraph::delete.edges(g, which(skip))
+        sp_verts <- suppressWarnings(igraph::shortest_paths(g1, s, t, weights = weights[!skip])$vpath[[1]])
+        if (length(sp_verts) < 2) {
+            skip[e] <- FALSE
+            next
+        }
+        sp_len <- path_length(sp_verts, xy)
+        if (sp_len >= max_distortion * elen[e]) {
+            skip[e] <- FALSE
+            next
+        }
+        lock[igraph::get.edge.ids(g, rep(as.integer(sp_verts), each = 2)[-c(1, 2 * length(sp_verts))])] <- TRUE
+        cpoints[[e]] <- xy[sp_verts, ]
     }
-    lock[igraph::get.edge.ids(g,rep(as.integer(sp_verts),each=2)[-c(1,2*length(sp_verts))])] <- TRUE
-    cpoints[[e]] <- xy[sp_verts,]
-  }
-  cpoints_bezier <- lapply(cpoints,approximateBezier,n=segments)
+    cpoints_bezier <- lapply(cpoints, approximateBezier, n = segments)
 
-  idx <- seq(0, 1, length.out = segments)
-  data_bundle <- as.data.frame(cbind(
-    do.call("rbind",cpoints_bezier),
-    rep(idx,m),
-    rep(1:m,each=segments)))
+    idx <- seq(0, 1, length.out = segments)
+    data_bundle <- as.data.frame(cbind(
+        do.call("rbind", cpoints_bezier),
+        rep(idx, m),
+        rep(1:m, each = segments)
+    ))
 
-  names(data_bundle) <- c("x","y","index","group")
+    names(data_bundle) <- c("x", "y", "index", "group")
 
-  data_bundle
+    data_bundle
 }
 
-path_length <- function(verts,xy){
-  plen <- 0
-  for(i in 1:(length(verts)-1)){
-    plen <- plen + sqrt((xy[i,1]-xy[i+1,1])^2+(xy[i,2]-xy[i+1,2])^2)
-  }
-  plen
+path_length <- function(verts, xy) {
+    plen <- 0
+    for (i in 1:(length(verts) - 1)) {
+        plen <- plen + sqrt((xy[i, 1] - xy[i + 1, 1])^2 + (xy[i, 2] - xy[i + 1, 2])^2)
+    }
+    plen
 }
 
 approximateBezier <- function(points, n) {
-  pnrow <- nrow(points)-1
-  tseq <- seq(0,1,length.out=n)
-  if(pnrow==1){
-    bezier <- cbind(
-      tseq*points[1,1]+(1-tseq)*points[2,1],
-      tseq*points[1,2]+(1-tseq)*points[2,2]
-      )
-  }
-  binoms <- choose(pnrow,seq(0,pnrow))
-  bezier <- matrix(0,length(tseq),2)
-  b <- 1
-  for(t in tseq){
-    p <- c(0,0)
-    for(i in 0:pnrow){
-      tpi <- (1-t)^(pnrow-i)
-      coeff <- tpi*t^i
-      p[1] <- p[1] + binoms[i+1] * coeff * points[i+1,1]
-      p[2] <- p[2] + binoms[i+1] * coeff * points[i+1,2]
+    pnrow <- nrow(points) - 1
+    tseq <- seq(0, 1, length.out = n)
+    if (pnrow == 1) {
+        bezier <- cbind(
+            tseq * points[1, 1] + (1 - tseq) * points[2, 1],
+            tseq * points[1, 2] + (1 - tseq) * points[2, 2]
+        )
+    }
+    binoms <- choose(pnrow, seq(0, pnrow))
+    bezier <- matrix(0, length(tseq), 2)
+    b <- 1
+    for (t in tseq) {
+        p <- c(0, 0)
+        for (i in 0:pnrow) {
+            tpi <- (1 - t)^(pnrow - i)
+            coeff <- tpi * t^i
+            p[1] <- p[1] + binoms[i + 1] * coeff * points[i + 1, 1]
+            p[2] <- p[2] + binoms[i + 1] * coeff * points[i + 1, 2]
+        }
+        bezier[b, ] <- p
+        b <- b + 1
     }
-    bezier[b,] <- p
-    b <- b+1
-  }
-  bezier
+    bezier
 }
diff --git a/R/bundle_stub.R b/R/bundle_stub.R
index f3405f1..09ec43a 100644
--- a/R/bundle_stub.R
+++ b/R/bundle_stub.R
@@ -15,235 +15,235 @@
 #' Nocaj, Arlind, and Ulrik Brandes. "Stub bundling and confluent spirals for geographic networks." International Symposium on Graph Drawing. Springer, Cham, 2013.
 #' @examples
 #' library(igraph)
-#' g <- graph.star(10,"undirected")
+#' g <- graph.star(10, "undirected")
 #'
 #' xy <- matrix(c(
-#'  0,0,
-#'  cos(90*pi/180),sin(90*pi/180),
-#'  cos(80*pi/180),sin(80*pi/180),
-#'  cos(70*pi/180),sin(70*pi/180),
-#'  cos(330*pi/180),sin(330*pi/180),
-#'  cos(320*pi/180),sin(320*pi/180),
-#'  cos(310*pi/180),sin(310*pi/180),
-#'  cos(210*pi/180),sin(210*pi/180),
-#'  cos(200*pi/180),sin(200*pi/180),
-#'  cos(190*pi/180),sin(190*pi/180)
-#'),ncol=2,byrow=TRUE)
+#'     0, 0,
+#'     cos(90 * pi / 180), sin(90 * pi / 180),
+#'     cos(80 * pi / 180), sin(80 * pi / 180),
+#'     cos(70 * pi / 180), sin(70 * pi / 180),
+#'     cos(330 * pi / 180), sin(330 * pi / 180),
+#'     cos(320 * pi / 180), sin(320 * pi / 180),
+#'     cos(310 * pi / 180), sin(310 * pi / 180),
+#'     cos(210 * pi / 180), sin(210 * pi / 180),
+#'     cos(200 * pi / 180), sin(200 * pi / 180),
+#'     cos(190 * pi / 180), sin(190 * pi / 180)
+#' ), ncol = 2, byrow = TRUE)
 #'
-#' edge_bundle_stub(g,xy)
+#' edge_bundle_stub(g, xy)
 #' # use ggforce::geom_bezier for plotting
 #' @export
 
-edge_bundle_stub <- function(object,xy,alpha = 11,beta = 75,gamma = 40,t = 0.5,tshift = 0.5){
-  if(any(class(object)=="igraph")){
-    if (!requireNamespace('igraph', quietly = TRUE)) {
-      stop('The `igraph` package is required for this functionality')
+edge_bundle_stub <- function(object, xy, alpha = 11, beta = 75, gamma = 40, t = 0.5, tshift = 0.5) {
+    if (any(class(object) == "igraph")) {
+        if (!requireNamespace("igraph", quietly = TRUE)) {
+            stop("The `igraph` package is required for this functionality")
+        }
+        el <- igraph::get.edgelist(object, FALSE)
+        adj <- igraph::get.adjlist(object, "all")
+    } else if (any(class(object) == "tbl_graph")) {
+        if (!requireNamespace("tidygraph", quietly = TRUE)) {
+            stop("The `tidygraph` package is required for this functionality")
+        }
+        object <- tidygraph::as.igraph(object)
+        el <- igraph::get.edgelist(object, FALSE)
+        adj <- igraph::get.adjlist(object, "all")
+    } else if (any(class(object) == "network")) {
+        el <- network::as.edgelist(object)
+
+        stop("`network` objects not supported. Convert object to `igraph` or `tbl_graph` first.")
+    } else {
+        stop("only `igraph` or `tbl_graph` objects supported.")
     }
-    el <-  igraph::get.edgelist(object,FALSE)
-    adj <- igraph::get.adjlist(object,"all")
 
-  } else if(any(class(object)=="tbl_graph")){
-    if (!requireNamespace('tidygraph', quietly = TRUE)) {
-      stop('The `tidygraph` package is required for this functionality')
+    idx <- seq(0, 1, length.out = 8) # 4 control points per stub/2 stubs
+
+    B <- compute_bundle_list(xy, adj, gamma, alpha)
+
+    res <- vector("list", nrow(el))
+    for (e in seq_len(nrow(el))) {
+        v <- el[e, 1]
+        w <- el[e, 2]
+        idw <- which(adj[[w]] == v)
+        Bw <- adj[[w]][which(B[[w]] == B[[w]][idw])]
+        idv <- which(adj[[v]] == w)
+        Bv <- adj[[v]][which(B[[v]] == B[[v]][idv])]
+        pv <- xy[v, ]
+        pw <- xy[w, ]
+        tst <- as.data.frame(control_points(xy, pv, pw, Bv, Bw, tshift, beta, t))
+        tst$idx <- idx
+        tst$grp <- c(rep(paste0(e, ".", 1), 4), rep(paste0(e, ".", 2), 4))
+        rownames(tst) <- NULL
+        res[[e]] <- tst
     }
-    object <- tidygraph::as.igraph(object)
-    el <-  igraph::get.edgelist(object,FALSE)
-    adj <- igraph::get.adjlist(object,"all")
-
-  } else if(any(class(object)=="network")){
-    el <- network::as.edgelist(object)
-
-    stop("`network` objects not supported. Convert object to `igraph` or `tbl_graph` first.")
-  } else{
-    stop("only `igraph` or `tbl_graph` objects supported.")
-  }
-
-  idx <- seq(0, 1, length.out = 8) #4 control points per stub/2 stubs
-
-  B <- compute_bundle_list(xy,adj,gamma,alpha)
-
-  res <- vector("list",nrow(el))
-  for(e in 1:nrow(el)){
-    v <- el[e,1]
-    w <- el[e,2]
-    idw <- which(adj[[w]]==v)
-    Bw <- adj[[w]][which(B[[w]]==B[[w]][idw])]
-    idv <- which(adj[[v]]==w)
-    Bv <- adj[[v]][which(B[[v]]==B[[v]][idv])]
-    pv <- xy[v,]
-    pw <- xy[w,]
-    tst <- as.data.frame(control_points(xy,pv,pw,Bv,Bw,tshift,beta,t))
-    tst$idx <- idx
-    tst$grp <- c(rep(paste0(e,".",1),4),rep(paste0(e,".",2),4))
-    rownames(tst) <- NULL
-    res[[e]] <- tst
-  }
-  data_bundle <- do.call(rbind,res)
-  names(data_bundle) <- c("x","y","index","group")
-
-  data_bundle
+    data_bundle <- do.call(rbind, res)
+    names(data_bundle) <- c("x", "y", "index", "group")
+
+    data_bundle
 }
 
 
-euclidean_dist <- function(p,q){
-  sqrt(sum((p-q)^2))
+euclidean_dist <- function(p, q) {
+    sqrt(sum((p - q)^2))
 }
 
-angle_edge <- function(P,Q){
-  pvec <- c(P[3]-P[1],P[4]-P[2])
-  qvec <- c(Q[3]-Q[1],Q[4]-Q[2])
-  # dot_pq <- sum(pvec*qvec)
-  # mag_pq <- sum(pvec^2)*sum(qvec^2)
-  # acos(dot_pq/mag_pq)*180/pi
-
-  m1 <- pvec[2]/pvec[1]
-  m2 <- qvec[2]/qvec[1]
-  a <- atan((m1 - m2 ) / (1+ m1*m2))
-  aa <- c(a*180/pi,-a*180/pi)
-  aa[aa<0] <- 360 + aa[aa<0]
-  aa <- aa*pi/180
-  min(aa)
+angle_edge <- function(P, Q) {
+    pvec <- c(P[3] - P[1], P[4] - P[2])
+    qvec <- c(Q[3] - Q[1], Q[4] - Q[2])
+    # dot_pq <- sum(pvec*qvec)
+    # mag_pq <- sum(pvec^2)*sum(qvec^2)
+    # acos(dot_pq/mag_pq)*180/pi
+
+    m1 <- pvec[2] / pvec[1]
+    m2 <- qvec[2] / qvec[1]
+    a <- atan((m1 - m2) / (1 + m1 * m2))
+    aa <- c(a * 180 / pi, -a * 180 / pi)
+    aa[aa < 0] <- 360 + aa[aa < 0]
+    aa <- aa * pi / 180
+    min(aa)
 }
 
-angle_edge_vec <- Vectorize( angle_edge )
-
-bundle_edges <- function(edges_xy,gamma,alpha){
-  if(nrow(edges_xy)==1){
-    return(1)
-    # return(t(c(1,1)))
-  }
-  dat <- as.data.frame(edges_xy)
-  dat[["V3"]] <- dat[["V3"]] - dat[["V1"]]
-  dat[["V4"]] <- dat[["V4"]] - dat[["V2"]]
-  dat[["V1"]] <- dat[["V2"]] <- 0
-  elen <- sqrt(dat[["V3"]]^2+dat[["V4"]]^2)
-  dat[["x"]] <- dat[["V3"]]/elen
-  dat[["y"]] <- dat[["V4"]]/elen
-  dat[["angle"]] <- atan2(dat[["y"]],dat[["x"]])*180/pi
-  dat[["angle"]] <- ifelse(dat[["angle"]]<0,360+dat[["angle"]],dat[["angle"]])
-  aord <- order(dat[["angle"]])
-  angle_vec <- dat[["angle"]][aord]
-  w <- pmin(abs(angle_vec-angle_vec[c(2:length(angle_vec),1)]),
-            360+angle_vec[c(2:length(angle_vec),1)]-angle_vec)
-
-  if(length(w)==1){
-    if(w>alpha){
-      return(c(1,2))
-    } else{
-      return(c(1,1))
+angle_edge_vec <- Vectorize(angle_edge)
+
+bundle_edges <- function(edges_xy, gamma, alpha) {
+    if (nrow(edges_xy) == 1) {
+        return(1)
+        # return(t(c(1,1)))
     }
-  }
-
-  start <- which.min(w)
-  bundles <- rep(0,length(aord))
-  if(start!=1){
-    w <- w[c(start:length(w),1:(start-1))]
-    aord <- aord[c(start:length(w),1:(start-1))]
-  }
-  bundles[aord[1]] <- 1
-  bundles[aord[2]] <- 1
-  cur <- 1
-  for(i in 2:length(w)-1){
-    if(w[i]<alpha & sum(bundles[bundles==cur])<gamma){
-      bundles[aord[i+1]] <- cur
-    } else{
-      cur <- cur+1
-      bundles[aord[i+1]] <- cur
+    dat <- as.data.frame(edges_xy)
+    dat[["V3"]] <- dat[["V3"]] - dat[["V1"]]
+    dat[["V4"]] <- dat[["V4"]] - dat[["V2"]]
+    dat[["V1"]] <- dat[["V2"]] <- 0
+    elen <- sqrt(dat[["V3"]]^2 + dat[["V4"]]^2)
+    dat[["x"]] <- dat[["V3"]] / elen
+    dat[["y"]] <- dat[["V4"]] / elen
+    dat[["angle"]] <- atan2(dat[["y"]], dat[["x"]]) * 180 / pi
+    dat[["angle"]] <- ifelse(dat[["angle"]] < 0, 360 + dat[["angle"]], dat[["angle"]])
+    aord <- order(dat[["angle"]])
+    angle_vec <- dat[["angle"]][aord]
+    w <- pmin(
+        abs(angle_vec - angle_vec[c(2:length(angle_vec), 1)]),
+        360 + angle_vec[c(2:length(angle_vec), 1)] - angle_vec
+    )
+
+    if (length(w) == 1) {
+        if (w > alpha) {
+            return(c(1, 2))
+        } else {
+            return(c(1, 1))
+        }
     }
-  }
 
-  # pa <- graph.ring(nrow(edges_xy),directed = FALSE)
-  # E(pa)$weight <- 360-w
-  # bundles <- cluster_louvain(pa,weights = E(pa)$weight)$membership[order(aord)]
+    start <- which.min(w)
+    bundles <- rep(0, length(aord))
+    if (start != 1) {
+        w <- w[c(start:length(w), 1:(start - 1))]
+        aord <- aord[c(start:length(w), 1:(start - 1))]
+    }
+    bundles[aord[1]] <- 1
+    bundles[aord[2]] <- 1
+    cur <- 1
+    for (i in 2:length(w) - 1) {
+        if (w[i] < alpha && sum(bundles[bundles == cur]) < gamma) {
+            bundles[aord[i + 1]] <- cur
+        } else {
+            cur <- cur + 1
+            bundles[aord[i + 1]] <- cur
+        }
+    }
 
-  bundles
+    # pa <- graph.ring(nrow(edges_xy),directed = FALSE)
+    # E(pa)$weight <- 360-w
+    # bundles <- cluster_louvain(pa,weights = E(pa)$weight)$membership[order(aord)]
+
+    bundles
 }
 
-compute_bundle_list <- function(xy,adj,gamma,alpha){
-  bundles_lst <-lapply(1:length(adj),function(i){
-    if(length(adj[[i]])>1){
-      edges_xy <- (cbind(matrix(xy[i,],nrow=length(adj[[i]]),ncol=2,byrow=TRUE),xy[adj[[i]],]))
-    }else{
-      edges_xy <- t((c(xy[i,],xy[adj[[i]],])))
-    }
-    bundle_edges(edges_xy,gamma,alpha)
-  })
-  bundles_lst
+compute_bundle_list <- function(xy, adj, gamma, alpha) {
+    bundles_lst <- lapply(seq_len(length(adj)), function(i) {
+        if (length(adj[[i]]) > 1) {
+            edges_xy <- (cbind(matrix(xy[i, ], nrow = length(adj[[i]]), ncol = 2, byrow = TRUE), xy[adj[[i]], ]))
+        } else {
+            edges_xy <- t((c(xy[i, ], xy[adj[[i]], ])))
+        }
+        bundle_edges(edges_xy, gamma, alpha)
+    })
+    bundles_lst
 }
 
-weighted_midpoint <- function(pv,pw,Bv,Bw,tshift){
-  0.5*(pv+pw)+(Bv/(Bv+Bw)-0.5)*tshift*(pw-pv)
+weighted_midpoint <- function(pv, pw, Bv, Bw, tshift) {
+    0.5 * (pv + pw) + (Bv / (Bv + Bw) - 0.5) * tshift * (pw - pv)
 }
 
-centroid <- function(Bv,xy){
-  if(length(Bv)>1){
-    colMeans(xy[Bv,])
-  } else{
-    xy[Bv,]
-  }
+centroid <- function(Bv, xy) {
+    if (length(Bv) > 1) {
+        colMeans(xy[Bv, ])
+    } else {
+        xy[Bv, ]
+    }
 }
 
-control_points <- function(xy,pv,pw,Bv,Bw,tshift,beta,t){
-  #pv2 is the point on c such that angle(pv,pv2,pm)=beta
-  beta <- beta*pi/180
+control_points <- function(xy, pv, pw, Bv, Bw, tshift, beta, t) {
+    # pv2 is the point on c such that angle(pv,pv2,pm)=beta
+    beta <- beta * pi / 180
 
-  #weighted midpoints
-  pm <- weighted_midpoint(pv,pw,length(Bv),length(Bw),tshift)
-  cv <- centroid(Bv,xy)
-  cw <- centroid(Bw,xy)
+    # weighted midpoints
+    pm <- weighted_midpoint(pv, pw, length(Bv), length(Bw), tshift)
+    cv <- centroid(Bv, xy)
+    cw <- centroid(Bw, xy)
 
-  #line between pv and cv
-  slope_v <- (cv[2] - pv[2])/(cv[1] - pv[1])
-  incep_v <-  cv[2] - slope_v * cv[1]
+    # line between pv and cv
+    slope_v <- (cv[2] - pv[2]) / (cv[1] - pv[1])
+    incep_v <- cv[2] - slope_v * cv[1]
 
-  #line between pw and cw
-  slope_w <- (cw[2] - pw[2])/(cw[1] - pw[1])
-  incep_w <-  cw[2] - slope_w * cw[1]
+    # line between pw and cw
+    slope_w <- (cw[2] - pw[2]) / (cw[1] - pw[1])
+    incep_w <- cw[2] - slope_w * cw[1]
 
-  # control points v -----------------------------------------------------------
-  #angle pm-pv-cv
-  alpha <- angle_edge(c(pv,pm),c(pv,cv))*pi/180
-  #remaining angle in triangle
-  gamma <- pi - alpha - beta
+    # control points v -----------------------------------------------------------
+    # angle pm-pv-cv
+    alpha <- angle_edge(c(pv, pm), c(pv, cv)) * pi / 180
+    # remaining angle in triangle
+    gamma <- pi - alpha - beta
 
-  #length of side opposite of beta
-  b <- euclidean_dist(pv,pm)
+    # length of side opposite of beta
+    b <- euclidean_dist(pv, pm)
 
-  l_pvcv <- sin(gamma) * b/sin(beta)
+    l_pvcv <- sin(gamma) * b / sin(beta)
 
-  pv2x <- c(pv[1] + l_pvcv/(sqrt(1+slope_v^2)),pv[1] - l_pvcv/(sqrt(1+slope_v^2)))
-  pv2y <- slope_v * pv2x + incep_v
+    pv2x <- c(pv[1] + l_pvcv / (sqrt(1 + slope_v^2)), pv[1] - l_pvcv / (sqrt(1 + slope_v^2)))
+    pv2y <- slope_v * pv2x + incep_v
 
-  idx <- which.min(c(euclidean_dist(pm,c(pv2x[1],pv2y[1])),euclidean_dist(pm,c(pv2x[2],pv2y[2]))))
-  pv2 <- c(pv2x[idx],pv2y[idx])
+    idx <- which.min(c(euclidean_dist(pm, c(pv2x[1], pv2y[1])), euclidean_dist(pm, c(pv2x[2], pv2y[2]))))
+    pv2 <- c(pv2x[idx], pv2y[idx])
 
-  pv1 <- pv + t * (pv2 - pv)
+    pv1 <- pv + t * (pv2 - pv)
 
-  # control points w -----------------------------------------------------------
-  #angle pm-pv-cv
-  alpha <- angle_edge(c(pw,pm),c(pw,cw))*pi/180
-  #remaining angle in triangle
-  gamma <- pi - alpha - beta
+    # control points w -----------------------------------------------------------
+    # angle pm-pv-cv
+    alpha <- angle_edge(c(pw, pm), c(pw, cw)) * pi / 180
+    # remaining angle in triangle
+    gamma <- pi - alpha - beta
 
-  #length of side opposite of beta
-  b <- euclidean_dist(pw,pm)
+    # length of side opposite of beta
+    b <- euclidean_dist(pw, pm)
 
-  l_pwcw <- sin(gamma) * b/sin(beta)
+    l_pwcw <- sin(gamma) * b / sin(beta)
 
-  pw2x <- c(pw[1] + l_pwcw/(sqrt(1+slope_w^2)),pw[1] - l_pwcw/(sqrt(1+slope_w^2)))
-  pw2y <- slope_w * pw2x + incep_w
+    pw2x <- c(pw[1] + l_pwcw / (sqrt(1 + slope_w^2)), pw[1] - l_pwcw / (sqrt(1 + slope_w^2)))
+    pw2y <- slope_w * pw2x + incep_w
 
-  idx <- which.min(c(euclidean_dist(pm,c(pw2x[1],pw2y[1])),euclidean_dist(pm,c(pw2x[2],pw2y[2]))))
-  pw2 <- c(pw2x[idx],pw2y[idx])
+    idx <- which.min(c(euclidean_dist(pm, c(pw2x[1], pw2y[1])), euclidean_dist(pm, c(pw2x[2], pw2y[2]))))
+    pw2 <- c(pw2x[idx], pw2y[idx])
 
-  pw1 <- pw + t * (pw2 - pw)
+    pw1 <- pw + t * (pw2 - pw)
 
-  # control point both
-  pvwm <- 0.5*(pv2 + pw2)
+    # control point both
+    pvwm <- 0.5 * (pv2 + pw2)
 
-  rbind(
-    rbind(pv,pv1,pv2,pvwm),
-    # rbind(pvwm,pw2,pw1,pw)
-    rbind(pw,pw1,pw2,pvwm)
-  )
+    rbind(
+        rbind(pv, pv1, pv2, pvwm),
+        # rbind(pvwm,pw2,pw1,pw)
+        rbind(pw, pw1, pw2, pvwm)
+    )
 }
diff --git a/R/convert_edges.R b/R/convert_edges.R
index 4afe78f..71335b0 100644
--- a/R/convert_edges.R
+++ b/R/convert_edges.R
@@ -6,54 +6,54 @@
 #' @author David Schoch
 #' @export
 #'
-convert_edges <- function(object,coords) UseMethod("convert_edges")
+convert_edges <- function(object, coords) UseMethod("convert_edges")
 
 #' @rdname convert_edges
 #' @method convert_edges default
 #' @export
-convert_edges.default <- function(object,coords){
-  stop("don't know how to handle class ", dQuote(data.class(object)))
+convert_edges.default <- function(object, coords) {
+    stop("don't know how to handle class ", dQuote(data.class(object)))
 }
 
 #' @rdname convert_edges
 #' @method convert_edges igraph
 #' @export
-convert_edges.igraph <- function(object,coords){
-  if (!requireNamespace('igraph', quietly = TRUE)) {
-    stop('The `igraph` package is required for this functionality')
-  }
-  el <- igraph::as_edgelist(object,names = FALSE)
-  if(igraph::vcount(object)!=nrow(coords)){
-    stop('number of rows in `coords` does not match number of vertices')
-  }
-  edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2])
-  edges_xy
+convert_edges.igraph <- function(object, coords) {
+    if (!requireNamespace("igraph", quietly = TRUE)) {
+        stop("The `igraph` package is required for this functionality")
+    }
+    el <- igraph::as_edgelist(object, names = FALSE)
+    if (igraph::vcount(object) != nrow(coords)) {
+        stop("number of rows in `coords` does not match number of vertices")
+    }
+    edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2])
+    edges_xy
 }
 
 #' @rdname convert_edges
 #' @method convert_edges network
 #' @export
-convert_edges.network <- function(object,coords){
-  if (!requireNamespace('network', quietly = TRUE)) {
-    stop('The `network` package is required for this functionality')
-  }
-  el <- network::as.edgelist(object)
-  if(network::get.network.attribute(object,"n")!=nrow(coords)){
-    stop('number of rows in `coords` does not match number of vertices')
-  }
-  edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2])
-  edges_xy
+convert_edges.network <- function(object, coords) {
+    if (!requireNamespace("network", quietly = TRUE)) {
+        stop("The `network` package is required for this functionality")
+    }
+    el <- network::as.edgelist(object)
+    if (network::get.network.attribute(object, "n") != nrow(coords)) {
+        stop("number of rows in `coords` does not match number of vertices")
+    }
+    edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2])
+    edges_xy
 }
 
 
 #' @rdname convert_edges
 #' @method convert_edges tbl_graph
 #' @export
-convert_edges.tbl_graph <- function(object,coords){
-  if (!requireNamespace('tidygraph', quietly = TRUE)) {
-    stop('The `tidygraph` package is required for this functionality')
-  }
-  el <- as.matrix(tidygraph::as_tibble(object,"edges"))
-  edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2])
-  edges_xy
+convert_edges.tbl_graph <- function(object, coords) {
+    if (!requireNamespace("tidygraph", quietly = TRUE)) {
+        stop("The `tidygraph` package is required for this functionality")
+    }
+    el <- as.matrix(tidygraph::as_tibble(object, "edges"))
+    edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2])
+    edges_xy
 }
diff --git a/R/data.R b/R/data.R
index 9e60e67..022fcaa 100644
--- a/R/data.R
+++ b/R/data.R
@@ -26,20 +26,3 @@
 #' @references
 #' Kujala, Rainer, et al. "A collection of public transport network data sets for 25 cities." Scientific data 5 (2018): 180089.
 "metro_berlin"
-
-# fl <- list.files("~/Documents/data/migration/",full.names = TRUE,pattern = "xls")
-# map(fl,function(f){
-#   df <- readxl::read_xls(f)
-#   names(df)[1] <- "first"
-#   idx <- which(df$first%in%state.name)
-#   idy <- which(df[6,]%in%state.name)
-#   tbl <- bind_cols(state.name,df[idx,idy])
-#   names(tbl) <- c("from",state.name)
-#   tbl <- tbl %>%
-#     gather("to","weight",Alabama:Wyoming) %>%
-#     mutate(weight=as.numeric(weight)) %>%
-#     dplyr::filter(!is.na(weight))  %>%
-#     mutate(year=parse_number(f))
-# }) -> tbl_lst
-#
-# us_migration <- as.data.frame(do.call(rbind,tbl_lst))
diff --git a/R/flow_tnss.R b/R/flow_tnss.R
index 2e702c7..e473a4c 100644
--- a/R/flow_tnss.R
+++ b/R/flow_tnss.R
@@ -33,54 +33,54 @@ tnss_dummies <- function(xy, root,
                          ndiag = 50,
                          ngrid = 50,
                          nrand = 50) {
-  n <- nrow(xy)
-  verts <- 1:n
-  leafs <- setdiff(verts, root)
-  dat <- matrix(0, 0, 2)
-
-  # circular points around leafs
-  if (circ) {
-    angles <- seq(0.01, 0.99 * 2 * pi, length.out = ncirc)
-    r <- rcirc
-    xy_circle <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] + r * cos(angles), xy[x, 2] + r * sin(angles))))
-    dat <- rbind(dat, xy_circle)
-  }
-
-  # points on line from source to leafs
-  if (line) {
-    tseq <- seq(0.2, 0.9, length.out = nline)
-    xy_lines <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] * tseq + xy[root, 1] * (1 - tseq), xy[x, 2] * tseq + xy[root, 2] * (1 - tseq))))
-    dat <- rbind(dat, xy_lines)
-  }
-
-  # diagonals through space
-  if (diag) {
-    pts_tr <- c(max(xy[, 1]), max(xy[, 2]))
-    pts_br <- c(max(xy[, 1]), min(xy[, 2]))
-    pts_bl <- c(min(xy[, 1]), min(xy[, 2]))
-    pts_tl <- c(max(xy[, 1]), max(xy[, 2]))
-    pts_extra <- rbind(pts_tr, pts_br, pts_bl, pts_tl)
-    tseq <- seq(0.1, 0.9, length.out = ndiag)
-    xy_extra <- do.call(rbind, lapply(1:4, function(x) cbind(pts_extra[x, 1] * tseq + xy[root, 1] * (1 - tseq), pts_extra[x, 2] * tseq + xy[root, 2] * (1 - tseq))))
-    dat <- rbind(dat, xy_extra)
-  }
-
-  # create an equidistant grid
-  if (grid) {
-    xdiff <- seq(min(xy[, 1]), max(xy[, 1]), length.out = ngrid)
-    ydiff <- seq(min(xy[, 2]), max(xy[, 2]), length.out = ngrid)
-    xy_grid <- as.matrix(expand.grid(xdiff, ydiff))
-    colnames(xy_grid) <- NULL
-    dat <- rbind(dat, xy_grid)
-  }
-
-  # some random points
-  if (rand) {
-    xy_rand <- cbind(stats::runif(nrand, min(xy[, 1]), max(xy[, 1])), stats::runif(50, min(xy[, 2]), max(xy[, 2])))
-    dat <- rbind(dat, xy_rand)
-  }
-
-  dat[!duplicated(dat), ]
+    n <- nrow(xy)
+    verts <- 1:n
+    leafs <- setdiff(verts, root)
+    dat <- matrix(0, 0, 2)
+
+    # circular points around leafs
+    if (circ) {
+        angles <- seq(0.01, 0.99 * 2 * pi, length.out = ncirc)
+        r <- rcirc
+        xy_circle <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] + r * cos(angles), xy[x, 2] + r * sin(angles))))
+        dat <- rbind(dat, xy_circle)
+    }
+
+    # points on line from source to leafs
+    if (line) {
+        tseq <- seq(0.2, 0.9, length.out = nline)
+        xy_lines <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] * tseq + xy[root, 1] * (1 - tseq), xy[x, 2] * tseq + xy[root, 2] * (1 - tseq))))
+        dat <- rbind(dat, xy_lines)
+    }
+
+    # diagonals through space
+    if (diag) {
+        pts_tr <- c(max(xy[, 1]), max(xy[, 2]))
+        pts_br <- c(max(xy[, 1]), min(xy[, 2]))
+        pts_bl <- c(min(xy[, 1]), min(xy[, 2]))
+        pts_tl <- c(max(xy[, 1]), max(xy[, 2]))
+        pts_extra <- rbind(pts_tr, pts_br, pts_bl, pts_tl)
+        tseq <- seq(0.1, 0.9, length.out = ndiag)
+        xy_extra <- do.call(rbind, lapply(1:4, function(x) cbind(pts_extra[x, 1] * tseq + xy[root, 1] * (1 - tseq), pts_extra[x, 2] * tseq + xy[root, 2] * (1 - tseq))))
+        dat <- rbind(dat, xy_extra)
+    }
+
+    # create an equidistant grid
+    if (grid) {
+        xdiff <- seq(min(xy[, 1]), max(xy[, 1]), length.out = ngrid)
+        ydiff <- seq(min(xy[, 2]), max(xy[, 2]), length.out = ngrid)
+        xy_grid <- as.matrix(expand.grid(xdiff, ydiff))
+        colnames(xy_grid) <- NULL
+        dat <- rbind(dat, xy_grid)
+    }
+
+    # some random points
+    if (rand) {
+        xy_rand <- cbind(stats::runif(nrand, min(xy[, 1]), max(xy[, 1])), stats::runif(50, min(xy[, 2]), max(xy[, 2])))
+        dat <- rbind(dat, xy_rand)
+    }
+
+    dat[!duplicated(dat), ]
 }
 
 #' @title Create Steiner tree from real and dummy points
@@ -104,121 +104,121 @@ tnss_dummies <- function(xy, root,
 #' @export
 
 tnss_tree <- function(g, xy, xydummy, root, gamma = 0.9, epsilon = 0.3, elen = Inf, order = "random") {
-  xymesh <- rbind(xy, xydummy)
-
-  n <- nrow(xy)
-  verts <- 1:n
-  leafs <- setdiff(verts, root)
-
-  # triangulate points
-  tria <- interp::tri.mesh(xymesh[, 1], xymesh[, 2], duplicate = "remove")
-
-  # create network
-  g1 <- igraph::graph_from_edgelist(rbind(tria$trlist[, 1:2], tria$trlist[, 2:3]), F)
-  g1 <- igraph::simplify(g1)
-  igraph::V(g1)$tnss <- "dummy"
-  igraph::V(g1)$tnss[seq_len(nrow(xy))] <- "real"
-  igraph::V(g1)$x <- xymesh[, 1]
-  igraph::V(g1)$y <- xymesh[, 2]
-
-  # delete edges that are too long
-  el <- igraph::get.edgelist(g1, FALSE)
-  edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2])
-  dist <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2))
-  idx <- which(dist > elen)
-  if (length(idx) != 0) {
-    g1 <- igraph::delete.edges(g1, idx)
-  }
-
-  # edge weights from distances
-  el <- igraph::get.edgelist(g1, FALSE)
-  edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2])
-  igraph::E(g1)$weight <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2))
-
-  # calculate all shortest paths to eliminate dummy nodes
-  sp_nodes <- vector("list", length(leafs))
-  sp_edges <- vector("list", length(leafs))
-  k <- 0
-
-  g2 <- igraph::as.directed(g1, "mutual")
-  ide <- which(igraph::get.edgelist(g2, FALSE)[, 1] %in% leafs)
-  g2 <- igraph::delete.edges(g2, ide)
-
-
-  dist2root <- sqrt((xy[root, 1] - xy[leafs, 1])^2 + (xy[root, 2] - xy[leafs, 2])^2)
-  # minew <- min(igraph::E(g2)$weight)
-  # leafs_order <- leafs[order(dist2root,decreasing = TRUE)]
-  if (order == "near") {
-    leafs <- leafs[order(dist2root, decreasing = FALSE)]
-  } else if (order == "far") {
-    leafs <- leafs[order(dist2root, decreasing = TRUE)]
-  } else if (order == "weight") {
-    leafs <- leafs
-  } else {
-    leafs <- sample(leafs)
-  }
-  for (dest in leafs) {
-    k <- k + 1
-    sp_list <- igraph::shortest_paths(g2, from = root, to = dest, weights = igraph::E(g2)$weight, output = "both")
-    sp_nodes[[k]] <- unlist(sp_list$vpath[[1]])
-    sp_edges[[k]] <- unlist(sp_list$epath[[1]])
-
-    igraph::E(g2)$weight[sp_edges[[k]]] <- gamma * igraph::E(g2)$weight[sp_edges[[k]]] #+0.01*minew
-  }
-  del_nodes <- unique(unlist(sp_nodes))
-  del_edges <- unique(unlist(sp_edges))
-
-  g3 <- igraph::delete.edges(g2, which(!((1:igraph::ecount(g2)) %in% del_edges)))
-  idx <- which(!igraph::V(g3) %in% del_nodes)
-  g3 <- igraph::delete.vertices(g3, idx)
-
-  # straighten edges
-  xymesh1 <- xymesh[-idx, ]
-  g4 <- igraph::as.undirected(g3)
-  g4 <- igraph::delete_edge_attr(g4, "weight")
-  deg <- igraph::degree(g4)
-  del2 <- c()
-  for (dest in leafs) {
-    sp <- unlist(igraph::shortest_paths(g4, from = root, to = dest)$vpath)
-    keep <- which(duplicated(
-      rbind(xymesh1[sp, ], visvalingam(xymesh1[sp, ], epsilon = epsilon)),
-      fromLast = TRUE
-    ))
-    del <- sp[-keep]
-    del <- del[deg[del] == 2]
-    del2 <- c(del2, del)
-  }
-
-  if (is.null(igraph::V(g4)$name)) {
-    igraph::V(g4)$name <- paste0("dummy_", 1:igraph::vcount(g4))
-  }
-  if (!is.null(igraph::V(g)$name)) {
-    igraph::V(g4)$name[1:n] <- igraph::V(g)$name
-  }
-  del2_name <- igraph::V(g4)$name[unique(del2)]
-  g5 <- g4
-  for (v in del2_name) {
-    ni <- igraph::neighborhood(g5, 1, v, mindist = 1)
-    g5 <- igraph::add.edges(g5, unlist(igraph::neighborhood(g5, 1, v, mindist = 1)))
-    g5 <- igraph::delete_vertices(g5, v)
-  }
-
-  # calculate flow from edge weight
-  gfinal <- g5
-  igraph::E(gfinal)$flow <- 0
-  el <- igraph::get.edgelist(g, FALSE)
-  for (dest in leafs) {
-    ide <- which(el[, 1] == dest | el[, 2] == dest)
-    w <- igraph::E(g)$weight[ide]
-    sp <- igraph::shortest_paths(gfinal, root, dest, weights = NA, output = "epath")
-
-    igraph::E(gfinal)$flow[unlist(sp$epath)] <- igraph::E(gfinal)$flow[unlist(sp$epath)] + w
-  }
-  gfinal$name <- "approx steiner tree"
-  igraph::V(gfinal)$tnss[root] <- "root"
-  igraph::V(gfinal)$tnss[leafs] <- "leaf"
-  class(gfinal) <- c("steiner_tree", class(gfinal))
-  gfinal
+    xymesh <- rbind(xy, xydummy)
+
+    n <- nrow(xy)
+    verts <- 1:n
+    leafs <- setdiff(verts, root)
+
+    # triangulate points
+    tria <- interp::tri.mesh(xymesh[, 1], xymesh[, 2], duplicate = "remove")
+
+    # create network
+    g1 <- igraph::graph_from_edgelist(rbind(tria$trlist[, 1:2], tria$trlist[, 2:3]), F)
+    g1 <- igraph::simplify(g1)
+    igraph::V(g1)$tnss <- "dummy"
+    igraph::V(g1)$tnss[seq_len(nrow(xy))] <- "real"
+    igraph::V(g1)$x <- xymesh[, 1]
+    igraph::V(g1)$y <- xymesh[, 2]
+
+    # delete edges that are too long
+    el <- igraph::get.edgelist(g1, FALSE)
+    edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2])
+    dist <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2))
+    idx <- which(dist > elen)
+    if (length(idx) != 0) {
+        g1 <- igraph::delete.edges(g1, idx)
+    }
+
+    # edge weights from distances
+    el <- igraph::get.edgelist(g1, FALSE)
+    edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2])
+    igraph::E(g1)$weight <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2))
+
+    # calculate all shortest paths to eliminate dummy nodes
+    sp_nodes <- vector("list", length(leafs))
+    sp_edges <- vector("list", length(leafs))
+    k <- 0
+
+    g2 <- igraph::as.directed(g1, "mutual")
+    ide <- which(igraph::get.edgelist(g2, FALSE)[, 1] %in% leafs)
+    g2 <- igraph::delete.edges(g2, ide)
+
+
+    dist2root <- sqrt((xy[root, 1] - xy[leafs, 1])^2 + (xy[root, 2] - xy[leafs, 2])^2)
+    # minew <- min(igraph::E(g2)$weight)
+    # leafs_order <- leafs[order(dist2root,decreasing = TRUE)]
+    if (order == "near") {
+        leafs <- leafs[order(dist2root, decreasing = FALSE)]
+    } else if (order == "far") {
+        leafs <- leafs[order(dist2root, decreasing = TRUE)]
+    } else if (order == "weight") {
+        leafs <- leafs
+    } else {
+        leafs <- sample(leafs)
+    }
+    for (dest in leafs) {
+        k <- k + 1
+        sp_list <- igraph::shortest_paths(g2, from = root, to = dest, weights = igraph::E(g2)$weight, output = "both")
+        sp_nodes[[k]] <- unlist(sp_list$vpath[[1]])
+        sp_edges[[k]] <- unlist(sp_list$epath[[1]])
+
+        igraph::E(g2)$weight[sp_edges[[k]]] <- gamma * igraph::E(g2)$weight[sp_edges[[k]]] #+0.01*minew
+    }
+    del_nodes <- unique(unlist(sp_nodes))
+    del_edges <- unique(unlist(sp_edges))
+
+    g3 <- igraph::delete.edges(g2, which(!((1:igraph::ecount(g2)) %in% del_edges)))
+    idx <- which(!igraph::V(g3) %in% del_nodes)
+    g3 <- igraph::delete.vertices(g3, idx)
+
+    # straighten edges
+    xymesh1 <- xymesh[-idx, ]
+    g4 <- igraph::as.undirected(g3)
+    g4 <- igraph::delete_edge_attr(g4, "weight")
+    deg <- igraph::degree(g4)
+    del2 <- c()
+    for (dest in leafs) {
+        sp <- unlist(igraph::shortest_paths(g4, from = root, to = dest)$vpath)
+        keep <- which(duplicated(
+            rbind(xymesh1[sp, ], visvalingam(xymesh1[sp, ], epsilon = epsilon)),
+            fromLast = TRUE
+        ))
+        del <- sp[-keep]
+        del <- del[deg[del] == 2]
+        del2 <- c(del2, del)
+    }
+
+    if (is.null(igraph::V(g4)$name)) {
+        igraph::V(g4)$name <- paste0("dummy_", 1:igraph::vcount(g4))
+    }
+    if (!is.null(igraph::V(g)$name)) {
+        igraph::V(g4)$name[1:n] <- igraph::V(g)$name
+    }
+    del2_name <- igraph::V(g4)$name[unique(del2)]
+    g5 <- g4
+    for (v in del2_name) {
+        ni <- igraph::neighborhood(g5, 1, v, mindist = 1)
+        g5 <- igraph::add.edges(g5, unlist(igraph::neighborhood(g5, 1, v, mindist = 1)))
+        g5 <- igraph::delete_vertices(g5, v)
+    }
+
+    # calculate flow from edge weight
+    gfinal <- g5
+    igraph::E(gfinal)$flow <- 0
+    el <- igraph::get.edgelist(g, FALSE)
+    for (dest in leafs) {
+        ide <- which(el[, 1] == dest | el[, 2] == dest)
+        w <- igraph::E(g)$weight[ide]
+        sp <- igraph::shortest_paths(gfinal, root, dest, weights = NA, output = "epath")
+
+        igraph::E(gfinal)$flow[unlist(sp$epath)] <- igraph::E(gfinal)$flow[unlist(sp$epath)] + w
+    }
+    gfinal$name <- "approx steiner tree"
+    igraph::V(gfinal)$tnss[root] <- "root"
+    igraph::V(gfinal)$tnss[leafs] <- "leaf"
+    class(gfinal) <- c("steiner_tree", class(gfinal))
+    gfinal
 }
 
 #' @title Smooth a Steiner tree
@@ -237,141 +237,141 @@ tnss_tree <- function(g, xy, xydummy, root, gamma = 0.9, epsilon = 0.3, elen = I
 #' @export
 
 tnss_smooth <- function(g, bw = 3, n = 10) {
-  if (!"steiner_tree" %in% class(g)) {
-    stop("g must be a steiner tree created with tnss_tree().")
-  }
+    if (!"steiner_tree" %in% class(g)) {
+        stop("g must be a steiner tree created with tnss_tree().")
+    }
 
-  root <- which(igraph::V(g)$tnss == "root")
-  leafs <- which(igraph::V(g)$tnss == "leaf")
+    root <- which(igraph::V(g)$tnss == "root")
+    leafs <- which(igraph::V(g)$tnss == "leaf")
 
-  el <- igraph::get.edgelist(g, names = FALSE)
-  xy <- cbind(igraph::V(g)$x, igraph::V(g)$y)
+    el <- igraph::get.edgelist(g, names = FALSE)
+    xy <- cbind(igraph::V(g)$x, igraph::V(g)$y)
 
-  ord <- order(igraph::distances(g, root, leafs), decreasing = TRUE)
-  res <- matrix(0, 0, 4)
+    ord <- order(igraph::distances(g, root, leafs), decreasing = TRUE)
+    res <- matrix(0, 0, 4)
 
-  for (v in ord) {
-    dest <- leafs[v]
-    sp <- unlist(igraph::shortest_paths(g, root, dest, output = "epath")$epath[[1]])
+    for (v in ord) {
+        dest <- leafs[v]
+        sp <- unlist(igraph::shortest_paths(g, root, dest, output = "epath")$epath[[1]])
 
-    path <- el[sp, , drop = FALSE]
-    edges_xy <- cbind(xy[path[, 1], 1], xy[path[, 1], 2], xy[path[, 2], 1], xy[path[, 2], 2])
-    pathxy <- matrix(c(t(edges_xy)), ncol = 2, byrow = TRUE)
-    pathxy <- pathxy[!duplicated(pathxy), ]
-    flow <- igraph::E(g)$flow[sp]
-    flow <- c(flow, flow[length(flow)])
-    pathxy <- cbind(pathxy, flow)
+        path <- el[sp, , drop = FALSE]
+        edges_xy <- cbind(xy[path[, 1], 1], xy[path[, 1], 2], xy[path[, 2], 1], xy[path[, 2], 2])
+        pathxy <- matrix(c(t(edges_xy)), ncol = 2, byrow = TRUE)
+        pathxy <- pathxy[!duplicated(pathxy), ]
+        flow <- igraph::E(g)$flow[sp]
+        flow <- c(flow, flow[length(flow)])
+        pathxy <- cbind(pathxy, flow)
 
-    if (nrow(pathxy) == 2) {
-      pathsm <- pathxy
-    } else {
-      pathsm <- edge_ksmooth(pathxy, bandwidth = bw, n = n)
-    }
+        if (nrow(pathxy) == 2) {
+            pathsm <- pathxy
+        } else {
+            pathsm <- edge_ksmooth(pathxy, bandwidth = bw, n = n)
+        }
 
-    res <- rbind(res, cbind(pathsm, dest))
-  }
-  df <- as.data.frame(res, row.names = NA)
-  names(df) <- c("x", "y", "flow", "destination")
-  df
+        res <- rbind(res, cbind(pathsm, dest))
+    }
+    df <- as.data.frame(res, row.names = NA)
+    names(df) <- c("x", "y", "flow", "destination")
+    df
 }
 
 # helpers ----
 DouglasPeucker <- function(points, epsilon) {
-  dmax <- 0
-  index <- 0
-  end <- nrow(points)
-  ResultList <- numeric(0)
-  if (end < 3) {
-    return(ResultList <- rbind(ResultList, points))
-  }
-  for (i in 2:(end - 1)) {
-    d <- ShortestDistance(points[i, ], line = rbind(points[1, ], points[end, ]))
-    if (d > dmax) {
-      index <- i
-      dmax <- d
+    dmax <- 0
+    index <- 0
+    end <- nrow(points)
+    ResultList <- numeric(0)
+    if (end < 3) {
+        return(ResultList <- rbind(ResultList, points))
+    }
+    for (i in 2:(end - 1)) {
+        d <- ShortestDistance(points[i, ], line = rbind(points[1, ], points[end, ]))
+        if (d > dmax) {
+            index <- i
+            dmax <- d
+        }
     }
-  }
-  # if dmax is greater than epsilon recursively apply
-  if (dmax > epsilon) {
-    recResults1 <- DouglasPeucker(points[1:index, ], epsilon)
-    recResults2 <- DouglasPeucker(points[index:end, ], epsilon)
-    ResultList <- rbind(ResultList, recResults1, recResults2)
-  } else {
-    ResultList <- rbind(ResultList, points[1, ], points[end, ])
-  }
-  ResultList <- as.matrix(ResultList[!duplicated(ResultList), ])
-  colnames(ResultList) <- c("x", "p")
-  return(ResultList)
+    # if dmax is greater than epsilon recursively apply
+    if (dmax > epsilon) {
+        recResults1 <- DouglasPeucker(points[1:index, ], epsilon)
+        recResults2 <- DouglasPeucker(points[index:end, ], epsilon)
+        ResultList <- rbind(ResultList, recResults1, recResults2)
+    } else {
+        ResultList <- rbind(ResultList, points[1, ], points[end, ])
+    }
+    ResultList <- as.matrix(ResultList[!duplicated(ResultList), ])
+    colnames(ResultList) <- c("x", "p")
+    return(ResultList)
 }
 
 ShortestDistance <- function(p, line) {
-  x1 <- line[1, 1]
-  y1 <- line[1, 2]
-  x2 <- line[2, 1]
-  y2 <- line[2, 2]
-  x0 <- p[1]
-  y0 <- p[2]
-  d <- abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt((y2 - y1)^2 + (x2 - x1)^2)
-  return(as.numeric(d))
+    x1 <- line[1, 1]
+    y1 <- line[1, 2]
+    x2 <- line[2, 1]
+    y2 <- line[2, 2]
+    x0 <- p[1]
+    y0 <- p[2]
+    d <- abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt((y2 - y1)^2 + (x2 - x1)^2)
+    return(as.numeric(d))
 }
 
 point_distance <- function(x) {
-  d <- diff(x)
-  sqrt(d[, 1]^2 + d[, 2]^2)
+    d <- diff(x)
+    sqrt(d[, 1]^2 + d[, 2]^2)
 }
 
 seq_multiple <- function(start, end, n) {
-  f <- function(x, y, z) {
-    sq <- seq(from = x, to = y, length.out = z)
-    sq[-1]
-  }
-  sq_mult <- mapply(f, start, end, n, SIMPLIFY = FALSE)
-  c(start[1], do.call(c, sq_mult))
+    f <- function(x, y, z) {
+        sq <- seq(from = x, to = y, length.out = z)
+        sq[-1]
+    }
+    sq_mult <- mapply(f, start, end, n, SIMPLIFY = FALSE)
+    c(start[1], do.call(c, sq_mult))
 }
 
 densify <- function(x, n = 10L) {
-  stopifnot(is.matrix(x), ncol(x) == 3, nrow(x) > 1)
-  n_pts <- nrow(x)
-  n_vec <- rep(n + 1, n_pts - 1)
-  x_dense <- seq_multiple(start = x[1:(n_pts - 1), 1], end = x[2:n_pts, 1], n = n_vec)
-  y_dense <- seq_multiple(start = x[1:(n_pts - 1), 2], end = x[2:n_pts, 2], n = n_vec)
-  f_dense <- seq_multiple(start = x[1:(n_pts - 1), 3], end = x[2:n_pts, 3], n = n_vec)
-  # f_dense <- c(rep(x[-nrow(x),3],each=n))
-  # f_dense <- c(f_dense,f_dense[length(f_dense)])
-
-  cbind(x_dense, y_dense, f_dense)
+    stopifnot(is.matrix(x), ncol(x) == 3, nrow(x) > 1)
+    n_pts <- nrow(x)
+    n_vec <- rep(n + 1, n_pts - 1)
+    x_dense <- seq_multiple(start = x[1:(n_pts - 1), 1], end = x[2:n_pts, 1], n = n_vec)
+    y_dense <- seq_multiple(start = x[1:(n_pts - 1), 2], end = x[2:n_pts, 2], n = n_vec)
+    f_dense <- seq_multiple(start = x[1:(n_pts - 1), 3], end = x[2:n_pts, 3], n = n_vec)
+    # f_dense <- c(rep(x[-nrow(x),3],each=n))
+    # f_dense <- c(f_dense,f_dense[length(f_dense)])
+
+    cbind(x_dense, y_dense, f_dense)
 }
 
 edge_ksmooth <- function(x, smoothness = 1, bandwidth = 2, n = 10L) {
-  pad <- list(
-    start = rbind(x[1, ], 2 * x[1, ] - x[2, ]),
-    end = rbind(x[nrow(x), ], 2 * x[nrow(x), ] - x[nrow(x) - 1, ])
-  )
-  pad$start[2, 3] <- pad$start[1, 3]
-
-  pad$start <- densify(pad$start[2:1, ], n = n)
-  pad$start[nrow(pad$start), 3] <- pad$start[nrow(pad$start) - 1, 3]
-  pad$end <- densify(pad$end, n = n)
-
-  x_dense <- densify(x, n = n)
-
-  n_pts <- nrow(x_dense)
-  x_pad <- rbind(pad$start, x_dense, pad$end)
-  d <- c(0, cumsum(point_distance(x_pad)))
-  x_smooth <- stats::ksmooth(d, x_pad[, 1],
-    n.points = length(d),
-    kernel = "normal", bandwidth = bandwidth
-  )
-  y_smooth <- stats::ksmooth(d, x_pad[, 2],
-    n.points = length(d),
-    kernel = "normal", bandwidth = bandwidth
-  )
-  keep_rows <- (x_smooth$x >= d[nrow(pad$start) + 1]) &
-    (x_smooth$x <= d[nrow(pad$start) + n_pts])
-  x_new <- cbind(x_smooth$y, y_smooth$y, x_pad[, 3])[keep_rows, ]
-  x_new[1, ] <- pad$start[nrow(pad$start), ]
-  x_new[nrow(x_new), ] <- pad$end[1, ]
-  return(x_new)
+    pad <- list(
+        start = rbind(x[1, ], 2 * x[1, ] - x[2, ]),
+        end = rbind(x[nrow(x), ], 2 * x[nrow(x), ] - x[nrow(x) - 1, ])
+    )
+    pad$start[2, 3] <- pad$start[1, 3]
+
+    pad$start <- densify(pad$start[2:1, ], n = n)
+    pad$start[nrow(pad$start), 3] <- pad$start[nrow(pad$start) - 1, 3]
+    pad$end <- densify(pad$end, n = n)
+
+    x_dense <- densify(x, n = n)
+
+    n_pts <- nrow(x_dense)
+    x_pad <- rbind(pad$start, x_dense, pad$end)
+    d <- c(0, cumsum(point_distance(x_pad)))
+    x_smooth <- stats::ksmooth(d, x_pad[, 1],
+        n.points = length(d),
+        kernel = "normal", bandwidth = bandwidth
+    )
+    y_smooth <- stats::ksmooth(d, x_pad[, 2],
+        n.points = length(d),
+        kernel = "normal", bandwidth = bandwidth
+    )
+    keep_rows <- (x_smooth$x >= d[nrow(pad$start) + 1]) &
+        (x_smooth$x <= d[nrow(pad$start) + n_pts])
+    x_new <- cbind(x_smooth$y, y_smooth$y, x_pad[, 3])[keep_rows, ]
+    x_new[1, ] <- pad$start[nrow(pad$start), ]
+    x_new[nrow(x_new), ] <- pad$end[1, ]
+    return(x_new)
 }
 
 
@@ -379,45 +379,45 @@ edge_ksmooth <- function(x, smoothness = 1, bandwidth = 2, n = 10L) {
 # credit to @coolbutuseless
 #' @importFrom utils head tail
 visvalingam <- function(points, epsilon) {
-  x <- points[, 1]
-  y <- points[, 2]
-  n <- max(c(floor(length(x) * epsilon), 2))
-  # Sanity Check
-  stopifnot(length(x) == length(y))
-  stopifnot(n <= length(x) && n >= 2)
-  if (length(x) == 2) {
-    return(list(x = x, y = y))
-  }
-
-  # Remove points
-  for (i in seq_len(length(x) - n)) {
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    # Find areas
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    x1 <- head(x, -2)
-    x2 <- head(x[-1], -1)
-    x3 <- tail(x, -2)
-
-    y1 <- head(y, -2)
-    y2 <- head(y[-1], -1)
-    y3 <- tail(y, -2)
-
-    a0 <- x1 - x2
-    a1 <- x3 - x2
-
-    b0 <- y1 - y2
-    b1 <- y3 - y2
-
-    tri_areas <- abs(a0 * b1 - a1 * b0) # / 2
-
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    # Find minimim area triangle and remove its centre vertex from all pts
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    min_tri <- which.min(tri_areas)
-
-    x <- x[-(min_tri + 1)]
-    y <- y[-(min_tri + 1)]
-  }
-
-  cbind(x, y)
+    x <- points[, 1]
+    y <- points[, 2]
+    n <- max(c(floor(length(x) * epsilon), 2))
+    # Sanity Check
+    stopifnot(length(x) == length(y))
+    stopifnot(n <= length(x) && n >= 2)
+    if (length(x) == 2) {
+        return(list(x = x, y = y))
+    }
+
+    # Remove points
+    for (i in seq_len(length(x) - n)) {
+        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        # Find areas
+        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        x1 <- head(x, -2)
+        x2 <- head(x[-1], -1)
+        x3 <- tail(x, -2)
+
+        y1 <- head(y, -2)
+        y2 <- head(y[-1], -1)
+        y3 <- tail(y, -2)
+
+        a0 <- x1 - x2
+        a1 <- x3 - x2
+
+        b0 <- y1 - y2
+        b1 <- y3 - y2
+
+        tri_areas <- abs(a0 * b1 - a1 * b0) # / 2
+
+        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        # Find minimim area triangle and remove its centre vertex from all pts
+        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        min_tri <- which.min(tri_areas)
+
+        x <- x[-(min_tri + 1)]
+        y <- y[-(min_tri + 1)]
+    }
+
+    cbind(x, y)
 }
diff --git a/R/metro_multicriteria.R b/R/metro_multicriteria.R
index ae7fd71..f9ecc76 100644
--- a/R/metro_multicriteria.R
+++ b/R/metro_multicriteria.R
@@ -21,60 +21,58 @@
 #' # the algorithm has problems with parallel edges
 #' library(igraph)
 #' g <- simplify(metro_berlin)
-#' xy <- cbind(V(g)$lon,V(g)$lat)*100
+#' xy <- cbind(V(g)$lon, V(g)$lat) * 100
 #'
 #' # the algorithm is not very stable. try playing with the parameters
-#' xy_new <- metro_multicriteria(g,xy,l = 2,gr = 0.5,w = c(100,100,1,1,100),bsize = 35)
+#' xy_new <- metro_multicriteria(g, xy, l = 2, gr = 0.5, w = c(100, 100, 1, 1, 100), bsize = 35)
 #' @export
-metro_multicriteria <- function(object,xy,l = 2,gr = 0.0025,w = rep(1,5),bsize = 5){
-  n <- igraph::vcount(object)
-  # lg <- l*gr
-  adj <- as_adj_list1(object)
-  adj <- lapply(adj, function(x) x-1)
-  adj_deg2 <- adj[unlist(lapply(adj,length))==2]
-  el <- igraph::get.edgelist(object,FALSE) - 1
+metro_multicriteria <- function(object, xy, l = 2, gr = 0.0025, w = rep(1, 5), bsize = 5) {
+    adj <- as_adj_list1(object)
+    adj <- lapply(adj, function(x) x - 1)
+    adj_deg2 <- adj[unlist(lapply(adj, length)) == 2]
+    el <- igraph::get.edgelist(object, FALSE) - 1
 
-  xy <- snap_to_grid(xy,gr)
+    xy <- snap_to_grid(xy, gr)
 
-  bbox <- station_bbox(xy,bsize,gr)
+    bbox <- station_bbox(xy, bsize, gr)
 
-  xy_new <- layout_as_metro_iter(adj,el,adj_deg2,xy,bbox,l,gr,w,bsize)
-  xy_new
+    xy_new <- layout_as_metro_iter(adj, el, adj_deg2, xy, bbox, l, gr, w, bsize)
+    xy_new
 }
 
-#helper ----
-snap_to_grid <- function(xy,gr){
-  xmin <- min(xy[,1])
-  xmax <- max(xy[,1])
-  ymin <- min(xy[,2])
-  ymax <- max(xy[,2])
+# helper ----
+snap_to_grid <- function(xy, gr) {
+    xmin <- min(xy[, 1])
+    xmax <- max(xy[, 1])
+    ymin <- min(xy[, 2])
+    ymax <- max(xy[, 2])
 
-  deltax <- seq(xmin-4*gr,xmax+4*gr,by=gr)
-  deltay <- seq(ymin-4*gr,ymax+4*gr,by=gr)
+    deltax <- seq(xmin - 4 * gr, xmax + 4 * gr, by = gr)
+    deltay <- seq(ymin - 4 * gr, ymax + 4 * gr, by = gr)
 
-  xdiff <- outer(xy[,1],deltax,function(x,y) abs(x-y))
-  ydiff <- outer(xy[,2],deltay,function(x,y) abs(x-y))
+    xdiff <- outer(xy[, 1], deltax, function(x, y) abs(x - y))
+    ydiff <- outer(xy[, 2], deltay, function(x, y) abs(x - y))
 
-  xy_new <- cbind(deltax[apply(xdiff,1,which.min)],deltay[apply(ydiff,1,which.min)])
-  dups <- duplicated(xy_new)
-  while(any(dups)){
-    xy_new[which(dups),] <- xy_new[which(dups),] + c(sample(c(1,-1),1)*gr,sample(c(1,-1),1)*gr)
+    xy_new <- cbind(deltax[apply(xdiff, 1, which.min)], deltay[apply(ydiff, 1, which.min)])
     dups <- duplicated(xy_new)
-  }
-  xy_new
+    while (any(dups)) {
+        xy_new[which(dups), ] <- xy_new[which(dups), ] + c(sample(c(1, -1), 1) * gr, sample(c(1, -1), 1) * gr)
+        dups <- duplicated(xy_new)
+    }
+    xy_new
 }
 
-station_bbox <- function(xy,bsize,gr){
-  cbind(xy - bsize * gr,xy + bsize * gr)
+station_bbox <- function(xy, bsize, gr) {
+    cbind(xy - bsize * gr, xy + bsize * gr)
 }
 
-as_adj_list1 <- function(g){
-  n <- igraph::vcount(g)
-  lapply(1:n,function(i){
-    x <- g[[i]][[1]]
-    attr(x,"env") <- NULL
-    attr(x,"graph") <- NULL
-    class(x) <- NULL
-    x
-  })
+as_adj_list1 <- function(g) {
+    n <- igraph::vcount(g)
+    lapply(1:n, function(i) {
+        x <- g[[i]][[1]]
+        attr(x, "env") <- NULL
+        attr(x, "graph") <- NULL
+        class(x) <- NULL
+        x
+    })
 }
diff --git a/cran-comments.md b/cran-comments.md
index 50dc3af..3c4a871 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,11 +1,10 @@
-# Update to 0.4.1
+# Update to 0.4.2
 
-* Fixed broken URL in README
-* bug fixes and non user facing changes
+* code clean and maintenance
 
 ## Test environments
 
-* Ubuntu 20.04: R 4.2.2
+* Ubuntu 22.04: R 4.3.2
 * win-builder (devel and release)
 
 ## R CMD check results
diff --git a/data/cali2010.rda b/data/cali2010.rda
index bd75d45..58bf33d 100644
Binary files a/data/cali2010.rda and b/data/cali2010.rda differ
diff --git a/data/metro_berlin.rda b/data/metro_berlin.rda
index 4e25389..2e6396a 100644
Binary files a/data/metro_berlin.rda and b/data/metro_berlin.rda differ
diff --git a/data/us_flights.rda b/data/us_flights.rda
index d5cc1e8..caab268 100644
Binary files a/data/us_flights.rda and b/data/us_flights.rda differ
diff --git a/man/edge_bundle_force.Rd b/man/edge_bundle_force.Rd
index 31776b5..9c56743 100644
--- a/man/edge_bundle_force.Rd
+++ b/man/edge_bundle_force.Rd
@@ -54,9 +54,14 @@ see \href{https://github.com/schochastics/edgebundle}{online} for plotting tips
 }
 \examples{
 library(igraph)
-g  <- graph_from_edgelist(matrix(c(1,12,2,11,3,10,4,9,5,8,6,7),ncol = 2,byrow = TRUE),FALSE)
-xy <- cbind(c(rep(0,6),rep(1,6)),c(1:6,1:6))
-edge_bundle_force(g,xy)
+g <- graph_from_edgelist(
+    matrix(c(
+        1, 12, 2, 11, 3, 10,
+        4, 9, 5, 8, 6, 7
+    ), ncol = 2, byrow = TRUE), FALSE
+)
+xy <- cbind(c(rep(0, 6), rep(1, 6)), c(1:6, 1:6))
+edge_bundle_force(g, xy)
 }
 \references{
 Holten, Danny, and Jarke J. Van Wijk. "Force-Directed Edge Bundling for Graph Visualization." Computer Graphics Forum (Blackwell Publishing Ltd) 28, no. 3 (2009): 983-990.
diff --git a/man/edge_bundle_path.Rd b/man/edge_bundle_path.Rd
index 699ecb3..a9ebe73 100644
--- a/man/edge_bundle_path.Rd
+++ b/man/edge_bundle_path.Rd
@@ -30,9 +30,12 @@ see \href{https://github.com/schochastics/edgebundle}{online} for plotting tips
 }
 \examples{
 library(igraph)
-g  <- graph_from_edgelist(matrix(c(1,2,1,6,1,4,2,3,3,4,4,5,5,6),ncol = 2,byrow = TRUE),FALSE)
-xy <- cbind(c(0,10,25,40,50,50),c(0,15,25,15,0,-10))
-edge_bundle_path(g,xy)
+g <- graph_from_edgelist(matrix(c(
+    1, 2, 1, 6,
+    1, 4, 2, 3, 3, 4, 4, 5, 5, 6
+), ncol = 2, byrow = TRUE), FALSE)
+xy <- cbind(c(0, 10, 25, 40, 50, 50), c(0, 15, 25, 15, 0, -10))
+edge_bundle_path(g, xy)
 }
 \references{
 Wallinger, M., Archambault, D., Auber, D., Nollenburg, M., & Peltonen, J. (2021). Edge-Path Bundling: A Less Ambiguous Edge Bundling Approach. IEEE Transactions on Visualization and Computer Graphics.
diff --git a/man/edge_bundle_stub.Rd b/man/edge_bundle_stub.Rd
index 694495d..bfe7fab 100644
--- a/man/edge_bundle_stub.Rd
+++ b/man/edge_bundle_stub.Rd
@@ -40,22 +40,22 @@ see \href{https://github.com/schochastics/edgebundle}{online} for plotting tips
 }
 \examples{
 library(igraph)
-g <- graph.star(10,"undirected")
+g <- graph.star(10, "undirected")
 
 xy <- matrix(c(
- 0,0,
- cos(90*pi/180),sin(90*pi/180),
- cos(80*pi/180),sin(80*pi/180),
- cos(70*pi/180),sin(70*pi/180),
- cos(330*pi/180),sin(330*pi/180),
- cos(320*pi/180),sin(320*pi/180),
- cos(310*pi/180),sin(310*pi/180),
- cos(210*pi/180),sin(210*pi/180),
- cos(200*pi/180),sin(200*pi/180),
- cos(190*pi/180),sin(190*pi/180)
-),ncol=2,byrow=TRUE)
+    0, 0,
+    cos(90 * pi / 180), sin(90 * pi / 180),
+    cos(80 * pi / 180), sin(80 * pi / 180),
+    cos(70 * pi / 180), sin(70 * pi / 180),
+    cos(330 * pi / 180), sin(330 * pi / 180),
+    cos(320 * pi / 180), sin(320 * pi / 180),
+    cos(310 * pi / 180), sin(310 * pi / 180),
+    cos(210 * pi / 180), sin(210 * pi / 180),
+    cos(200 * pi / 180), sin(200 * pi / 180),
+    cos(190 * pi / 180), sin(190 * pi / 180)
+), ncol = 2, byrow = TRUE)
 
-edge_bundle_stub(g,xy)
+edge_bundle_stub(g, xy)
 # use ggforce::geom_bezier for plotting
 }
 \references{
diff --git a/man/edgebundle-package.Rd b/man/edgebundle-package.Rd
index e6c7a01..a6494b1 100644
--- a/man/edgebundle-package.Rd
+++ b/man/edgebundle-package.Rd
@@ -20,8 +20,8 @@ It includes:
 \seealso{
 Useful links:
 \itemize{
-  \item \url{http://edgebundle.schochastics.net/}
   \item \url{https://github.com/schochastics/edgebundle}
+  \item \url{https://schochastics.github.io/edgebundle/}
   \item Report bugs at \url{https://github.com/schochastics/edgebundle/issues}
 }
 
diff --git a/man/metro_multicriteria.Rd b/man/metro_multicriteria.Rd
index 40418c7..2a49531 100644
--- a/man/metro_multicriteria.Rd
+++ b/man/metro_multicriteria.Rd
@@ -40,10 +40,10 @@ see \href{https://github.com/schochastics/edgebundle}{online} for more plotting
 # the algorithm has problems with parallel edges
 library(igraph)
 g <- simplify(metro_berlin)
-xy <- cbind(V(g)$lon,V(g)$lat)*100
+xy <- cbind(V(g)$lon, V(g)$lat) * 100
 
 # the algorithm is not very stable. try playing with the parameters
-xy_new <- metro_multicriteria(g,xy,l = 2,gr = 0.5,w = c(100,100,1,1,100),bsize = 35)
+xy_new <- metro_multicriteria(g, xy, l = 2, gr = 0.5, w = c(100, 100, 1, 1, 100), bsize = 35)
 }
 \references{
 Stott, Jonathan, et al. "Automatic metro map layout using multicriteria optimization." IEEE Transactions on Visualization and Computer Graphics 17.1 (2010): 101-114.
diff --git a/man/tnss_dummies.Rd b/man/tnss_dummies.Rd
index 7dee6c9..64c8d9d 100644
--- a/man/tnss_dummies.Rd
+++ b/man/tnss_dummies.Rd
@@ -55,8 +55,8 @@ uses various sampling strategies to create dummy nodes for the \link{tnss_tree}
 }
 \examples{
 # dummy nodes for tree rooted in California
-xy <- cbind(state.center$x,state.center$y)
-xy_dummy <- tnss_dummies(xy,4)
+xy <- cbind(state.center$x, state.center$y)
+xy_dummy <- tnss_dummies(xy, 4)
 }
 \author{
 David Schoch
diff --git a/man/tnss_smooth.Rd b/man/tnss_smooth.Rd
index e8b3706..06b10fa 100644
--- a/man/tnss_smooth.Rd
+++ b/man/tnss_smooth.Rd
@@ -23,10 +23,10 @@ Converts the Steiner tree to smooth paths
 see see \href{https://github.com/schochastics/edgebundle}{online} for tips on plotting the result
 }
 \examples{
-xy <- cbind(state.center$x,state.center$y)[!state.name\%in\%c("Alaska","Hawaii"),]
-xy_dummy <- tnss_dummies(xy,root = 4)
-gtree <- tnss_tree(cali2010,xy,xy_dummy,root = 4,gamma = 0.9)
-tree_smooth <- tnss_smooth(gtree,bw = 10,n = 10)
+xy <- cbind(state.center$x, state.center$y)[!state.name \%in\% c("Alaska", "Hawaii"), ]
+xy_dummy <- tnss_dummies(xy, root = 4)
+gtree <- tnss_tree(cali2010, xy, xy_dummy, root = 4, gamma = 0.9)
+tree_smooth <- tnss_smooth(gtree, bw = 10, n = 10)
 }
 \author{
 David Schoch
diff --git a/man/tnss_tree.Rd b/man/tnss_tree.Rd
index 51599e8..4929bee 100644
--- a/man/tnss_tree.Rd
+++ b/man/tnss_tree.Rd
@@ -42,9 +42,9 @@ creates an approximated Steiner tree for a flow map visualization
 Use \link{tnss_smooth} to smooth the edges of the tree
 }
 \examples{
-xy <- cbind(state.center$x,state.center$y)[!state.name\%in\%c("Alaska","Hawaii"),]
-xy_dummy <- tnss_dummies(xy,root = 4)
-gtree <- tnss_tree(cali2010,xy,xy_dummy,root = 4,gamma = 0.9)
+xy <- cbind(state.center$x, state.center$y)[!state.name \%in\% c("Alaska", "Hawaii"), ]
+xy_dummy <- tnss_dummies(xy, root = 4)
+gtree <- tnss_tree(cali2010, xy, xy_dummy, root = 4, gamma = 0.9)
 }
 \references{
 Sun, Shipeng. "An automated spatial flow layout algorithm using triangulation, approximate Steiner tree, and path smoothing." AutoCarto, 2016.