diff --git a/DESCRIPTION b/DESCRIPTION
index 90ccbaf..40c78d1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: archeofrag
 Type: Package
 Title: Refitting and Spatial Analysis in Archaeology
-Version: 0.8
-Date: 2022-03-08
+Version: 0.8.1
+Date: 2022-07-13
 Authors@R: person("Sebastien", "Plutniak", email = "sebastien.plutniak@posteo.net", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6674-3806"))
 Author: Sebastien Plutniak [aut, cre] (<https://orcid.org/0000-0002-6674-3806>)
 Maintainer: Sebastien Plutniak <sebastien.plutniak@posteo.net>
@@ -10,14 +10,15 @@ Description: Methods to analyse fragmented objects in archaeology using refittin
 License: GPL-3
 Repository: CRAN
 Encoding: UTF-8
-Depends: igraph, RBGL
 Imports:
+	igraph,
     graphics,
     stats,
     grDevices,
     methods,
     utils
 Suggests: 
+	RBGL,
     knitr,
     rmarkdown,
     markdown
diff --git a/NAMESPACE b/NAMESPACE
index 8f8e219..b6c4c45 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,11 +18,52 @@ export(make_sr_graph)
 export(make_crsr_graph)
 export(make_frag_object)
 exportClasses(Frag.object)
-import(igraph)
-importFrom(grDevices, rainbow, rgb)
+importFrom(grDevices, rgb)
 importFrom(graphics, plot, hist)
 importFrom(methods, new)
 importFrom(utils, combn)
-importFrom(stats, na.omit, sd)
-importFrom(stats, wilcox.test)
-importFrom(RBGL, boyerMyrvoldPlanarityTest)
+importFrom(stats, na.omit, sd, wilcox.test)
+importFrom(igraph, 
+	add_edges, 
+	add_layout_, 
+	add_vertices, 
+	as_edgelist, 
+	as_graphnel, 
+	bibcoupling, 
+	components, 
+	component_wise, 
+	complementer,
+	"%--%",
+	decompose, 
+	degree, 
+	delete_vertex_attr, 
+	delete_vertices, 
+	diameter, 
+	distance_table, 
+	E, 
+	edge_attr, 
+	ego,
+	gorder, 
+	graph_attr, 
+	graph_from_adjacency_matrix, 
+	graph_from_data_frame, 
+	gsize, 
+	induced_subgraph, 
+	is_igraph, 
+	is_directed,
+	layout_with_fr, 
+	make_ring, 
+	make_graph,
+	plot.igraph, 
+	set_graph_attr, 
+	set_vertex_attr, 
+	simplify, 
+	subgraph.edges, 
+	subgraph_isomorphisms, 
+	transitivity, 
+	union,
+	V, 
+	vertex_attr, 
+	vertex_attr_names, 
+	with_fr) 
+
diff --git a/NEWS.md b/NEWS.md
index 15a1e31..313d22a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,11 @@
 
+# archeofrag 0.8.1
+Released: 2022-07-13
+
+* all functions from other packages are tagged with their package.
+* the RBGL is moved to the 'suggested' packages list, the functionalities related to planarity are now optional.
+* igraph functions' names are updated.
+* creation of a utils.R file, including internal functions to check the fragmentation graph and their "layer" argument.
 
 # archeofrag 0.8
 Released: 2022-03-08
diff --git a/R/frag.cycles.R b/R/frag.cycles.R
index 7d3f0af..82d8512 100644
--- a/R/frag.cycles.R
+++ b/R/frag.cycles.R
@@ -1,8 +1,8 @@
 
 
 .detect.cycle <- function(graph, x){
-  pattern <- make_ring(x) # create the pattern to detect:
-  cycles.list <- subgraph_isomorphisms(pattern, graph)
+  pattern <- igraph::make_ring(x) # create the pattern to detect:
+  cycles.list <- igraph::subgraph_isomorphisms(pattern, graph)
   res <- data.frame()
   if(length(cycles.list) != 0){ 
     res <- do.call(rbind, cycles.list) # convert the list to a data frame
@@ -19,10 +19,12 @@
   .max.cycles.only(x-1, res) # recursive call
 }
 
-#'	@export
+# 
 frag.cycles <- function(graph, kmax, max.cycles.only=FALSE) {
-  if (! is.igraph(graph)) stop("Not an igraph object")
-  if(graph_attr(graph, "frag_type") == "connection and similarity" ){
+  # tests:
+  .check.frag.graph(graph)  
+  
+  if(igraph::graph_attr(graph, "frag_type") == "connection and similarity" ){
     warning("Cycle detection in a 'connection and similarity' fragmentation graph is meaningless.")
   }
   if( ! is.numeric(kmax)){
diff --git a/R/frag.diameters.R b/R/frag.diameters.R
index d12c13c..9bccfad 100644
--- a/R/frag.diameters.R
+++ b/R/frag.diameters.R
@@ -1,18 +1,20 @@
 
-#'	@export
-frag.diameters <- function(graph, cumulative=FALSE) 
-{
-  if( ! is.igraph(graph)) stop("Not a graph object")
-  if( ! is.null(graph_attr(graph, "frag_type")) ) {
-    if(graph_attr(graph, "frag_type") == "connection and similarity"){
-      warning("Diameter distribution is meaningless for 'connection and similarity' fragmentation graphs")
+#
+frag.diameters <- function(graph, cumulative=FALSE){
+  # tests:
+  .check.frag.graph(graph)
+  
+  if( ! is.null(igraph::graph_attr(graph, "frag_type")) ) {
+    if(igraph::graph_attr(graph, "frag_type") == "connection and similarity"){
+      warning("Diameter distribution is meaningless for 'connection and similarity' fragmentation graphs.")
     }
   }
   if(! is.logical(cumulative)) stop("The 'cumulative' parameter must be logical.")
   
-  g.list <- decompose(graph)
-  diameter.v <- sapply(g.list, diameter, weights=NA, directed=FALSE, unconnected=FALSE)
-  diameter.v <- hist(diameter.v, breaks=0:max(diameter.v), plot=FALSE)$count 
+  # main body:
+  g.list <- igraph::decompose(graph)
+  diameter.v <- sapply(g.list, igraph::diameter, weights=NA, directed=FALSE, unconnected=FALSE)
+  diameter.v <- graphics::hist(diameter.v, breaks=0:max(diameter.v), plot=FALSE)$count 
   if(cumulative){
     diameter.v <- rev(cumsum(rev(diameter.v))) / sum(diameter.v)
   }
diff --git a/R/frag.edges.weighting.R b/R/frag.edges.weighting.R
index 1f0258f..e3ee05b 100644
--- a/R/frag.edges.weighting.R
+++ b/R/frag.edges.weighting.R
@@ -1,9 +1,9 @@
 .euclidean.distance <- function(graph, x, y, z){
-  coords <- cbind(name = V(graph)$name,
-                  x = get.vertex.attribute(graph, x),
-                  y = get.vertex.attribute(graph, y),
-                  z = get.vertex.attribute(graph, z))
-  e.list <- get.edgelist(graph)
+  coords <- cbind(name = igraph::V(graph)$name,
+                  x = igraph::vertex_attr(graph, x),
+                  y = igraph::vertex_attr(graph, y),
+                  z = igraph::vertex_attr(graph, z))
+  e.list <- igraph::as_edgelist(graph)
   colnames(e.list) <- c("name1", "name2")
   e.list <- merge(e.list, coords, by.x="name1", by.y="name", sort=FALSE)
   e.list <- merge(e.list, coords, by.x="name2", by.y="name",
@@ -27,12 +27,12 @@
 .get.morpho.spatial.params <- function(graph, x, y, z){
   # input: a graph and the names for x y z attributes
   # output: a data frame with the distance and morpho parameters for the graph edges
-  e.list <- get.edgelist(graph)
+  e.list <- igraph::as_edgelist(graph)
   e.list <- data.frame(e.list) 
   e.list$id <- 1:nrow(e.list) # a an it to sort the row later
   colnames(e.list) <- c("name1", "name2", "id")
-  morpho.df <- cbind(name = V(graph)$name, 
-                     morphometry = V(graph)$morphometry)
+  morpho.df <- cbind(name = igraph::V(graph)$name, 
+                     morphometry = igraph::V(graph)$morphometry)
   e.list <- merge(e.list, morpho.df, by.x="name1", by.y="name", sort=FALSE)
   e.list <- merge(e.list, morpho.df, by.x="name2", by.y="name",
                   suffixes=c("1", "2"), sort=FALSE)
@@ -49,7 +49,7 @@
   
   # compute spatial distances (this is done in this function to avoid repeating 3 times the command)
   e.list$distance <- 1 # default value
-  coords <- c(x, y, z)[ c(x, y, z) %in% vertex_attr_names(graph) ]
+  coords <- c(x, y, z)[ c(x, y, z) %in% igraph::vertex_attr_names(graph) ]
   if( length(coords) > 1) {   # if at least two coordinates
     e.list$distance <- .euclidean.distance(graph, x, y, z) + 1 # +1 to ensure a neutral value when all distances=0
   }
@@ -67,12 +67,12 @@
 
 .set.edge.weight <- function(g){
   #edges weight = sum of the degrees * transitivity-based index
-  if(gsize(g) == 0 )  return(0)
-  e.list <- as_edgelist(g)
+  if(igraph::gsize(g) == 0 )  return(0)
+  e.list <- igraph::as_edgelist(g)
   colnames(e.list) <- c("src", "tgt")
   v.list <- cbind(name = V(g)$name,
                   deg = igraph::degree(g),
-                  trans = transitivity(g, type = "localundirected",
+                  trans = igraph::transitivity(g, type = "localundirected",
                                        isolates = "zero", weights = NULL))
   e.list <- merge(e.list, v.list, by.x="src", by.y="name", sort=FALSE)
   e.list <- merge(e.list, v.list, by.x="tgt", by.y="name", sort=FALSE)
@@ -86,20 +86,11 @@
   e.list$SumDegree * e.list$trans.factor * size.factor
 }
 
-frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="", z="") 
-{
-  if(! is.igraph(graph)){
-    stop("Not a graph object")
-  }
-  if(! is.character(layer.attr)) {
-    stop("A character value is expected for the 'layer.attr' parameter.")
-  }
-  if (is.null(layer.attr)){
-    stop("No 'layer.attr' argument")
-  }
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
+frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="", z=""){
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
+  
   if(! is.character(morphometry)) {
     stop("A character value is expected for the 'morphometry' parameter.")
   }
@@ -107,28 +98,28 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
     stop("Character values are expected for the 'x', 'y' and 'z' parameters.")
   }
   
-  if (is.null(V(graph)$name)){
-    V(graph)$name <- seq(1:gorder(graph))
+  if(is.null(V(graph)$name)){
+    igraph::V(graph)$name <- seq(1:gorder(graph))
   }
   
-  missing.coords <- c(x, y, z)[ ! c(x, y, z) %in% c("", vertex_attr_names(graph)) ]
+  missing.coords <- c(x, y, z)[ ! c(x, y, z) %in% c("", igraph::vertex_attr_names(graph)) ]
   if(length(missing.coords) > 0){
     warning(paste("Missing coordinates:", paste(missing.coords, collapse = ", ")  ))
   }
   
-  V(graph)$morphometry <- 1 # set default value
+  igraph::V(graph)$morphometry <- 1 # set default value
   if( morphometry != ""){
-    if( morphometry %in% names(vertex_attr(graph)) ){
-      V(graph)$morphometry <- vertex_attr(graph, morphometry)
+    if( morphometry %in% names(igraph::vertex_attr(graph)) ){
+      igraph::V(graph)$morphometry <- igraph::vertex_attr(graph, morphometry)
     }else{
       stop(paste("No '", morphometry, "' vertices attribute", sep=""))
     }
   }
-  if(! is.numeric(V(graph)$morphometry)){
+  if(! is.numeric(igraph::V(graph)$morphometry)){
     stop("Numeric values are required for the 'morphometry' parameter.")
   }
   
-  layers <- vertex_attr(graph, layer.attr)
+  layers <- igraph::vertex_attr(graph, layer.attr)
   layers.u <- unique(layers)
   if(length(layers.u) > 2){
     stop("There are more than two layers.")
@@ -136,7 +127,7 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
   
   v1 <- layers == layers.u[1]
   v2 <- layers == layers.u[2]
-  E(graph)$id <- 1:gsize(graph)
+  igraph::E(graph)$id <- 1:igraph::gsize(graph)
   
   # get the max morphometric/spatial values observed in the data set:
   params  <- .get.morpho.spatial.params(graph, x, y, z)
@@ -147,23 +138,23 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
   
   # add a "distance" edge attribute if necessary:
   if(! is.null(params$distance)){
-    E(graph)$distance <- params$distance
+    igraph::E(graph)$distance <- params$distance
   }
   
   # generate subgraphs
-  g1  <- subgraph.edges(graph, E(graph)[ V(graph)[v1] %--% V(graph)[v1] ])
-  g2  <- subgraph.edges(graph, E(graph)[ V(graph)[v2] %--% V(graph)[v2] ])
-  g12 <- subgraph.edges(graph, E(graph)[ V(graph)[v1] %--% V(graph)[v2] ])
+  g1  <- igraph::subgraph.edges(graph, igraph::E(graph)[ igraph::V(graph)[v1] %--% igraph::V(graph)[v1] ])
+  g2  <- igraph::subgraph.edges(graph, igraph::E(graph)[ igraph::V(graph)[v2] %--% igraph::V(graph)[v2] ])
+  g12 <- igraph::subgraph.edges(graph, igraph::E(graph)[ igraph::V(graph)[v1] %--% igraph::V(graph)[v2] ])
   
   # extract edges indices:
-  e1  <- E(graph)$id %in% E(g1)$id
-  e2  <- E(graph)$id %in% E(g2)$id
-  e12 <- E(graph)$id %in% E(g12)$id
+  e1  <- igraph::E(graph)$id %in% igraph::E(g1)$id
+  e2  <- igraph::E(graph)$id %in% igraph::E(g2)$id
+  e12 <- igraph::E(graph)$id %in% igraph::E(g12)$id
   
   # get the morphometric/spatial parameters of the edges compute the 
   # morphometric/spatial factor for the edges  of each subgraphs:
   morpho.spatial.factor.g1 <- 0
-  if(gsize(g1) > 0){
+  if(igraph::gsize(g1) > 0){
     params1  <- .get.morpho.spatial.params(g1, x, y, z)
     morpho.spatial.factor.g1 <- apply(params1, 1, function(x){
       .morpho.spatial.factor(morpho1 = x[1],
@@ -177,7 +168,7 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
   }
   
   morpho.spatial.factor.g2 <- 0
-  if(gsize(g2) > 0){
+  if(igraph::gsize(g2) > 0){
     params2  <- .get.morpho.spatial.params(g2, x, y, z)
     morpho.spatial.factor.g2 <- apply(params2, 1, function(x){
       .morpho.spatial.factor(morpho1 = x[1],
@@ -191,7 +182,7 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
   }
   
   morpho.spatial.factor.g12 <- 0
-  if(gsize(g12) > 0){
+  if(igraph::gsize(g12) > 0){
     params12 <- .get.morpho.spatial.params(g12, x, y, z)
     morpho.spatial.factor.g12 <- apply(params12, 1, function(x){
       .morpho.spatial.factor(morpho1 = x[1],
@@ -205,19 +196,19 @@ frag.edges.weighting <- function(graph, layer.attr, morphometry="", x="", y="",
   }
   
   # compute the weights for the edges of the 3 subsets:
-  E(graph)[ e1 ]$weight <- .set.edge.weight(g1)
-  E(graph)[ e2 ]$weight <- .set.edge.weight(g2)
-  E(graph)[ e12]$weight <- .set.edge.weight(g12)
+  igraph::E(graph)[ e1 ]$weight <- .set.edge.weight(g1)
+  igraph::E(graph)[ e2 ]$weight <- .set.edge.weight(g2)
+  igraph::E(graph)[ e12]$weight <- .set.edge.weight(g12)
   #   apply the morpho/spatial modifier is required:
   if(sum(sapply(list(morphometry, x, y, z), function(x) x != ""))){
-    E(graph)[ e1 ]$weight <- E(graph)[e1]$weight * morpho.spatial.factor.g1
-    E(graph)[ e2 ]$weight <- E(graph)[e2]$weight * morpho.spatial.factor.g2
-    E(graph)[ e12]$weight <- E(graph)[e12]$weight* morpho.spatial.factor.g12
+    igraph::E(graph)[ e1 ]$weight <- igraph::E(graph)[e1]$weight * morpho.spatial.factor.g1
+    igraph::E(graph)[ e2 ]$weight <- igraph::E(graph)[e2]$weight * morpho.spatial.factor.g2
+    igraph::E(graph)[ e12]$weight <- igraph::E(graph)[e12]$weight* morpho.spatial.factor.g12
   }  
   
   # add tags to the edges:
-  E(graph)$scope <- "intra"
-  E(graph)[e12]$scope <- "extra"
+  igraph::E(graph)$scope <- "intra"
+  igraph::E(graph)[e12]$scope <- "extra"
   
   graph
 }
diff --git a/R/frag.get.layers.R b/R/frag.get.layers.R
index fe526ab..bcaf302 100644
--- a/R/frag.get.layers.R
+++ b/R/frag.get.layers.R
@@ -1,17 +1,15 @@
 
 frag.get.layers <- function(graph, layer.attr, sel.layers){
-  if(! is.igraph(graph))  stop("Not a graph object")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
   
-  layers <- vertex_attr(graph, layer.attr)
+  # function:
+  layers <- igraph::vertex_attr(graph, layer.attr)
   if(sum(sel.layers %in% layers) != length(sel.layers)) stop("Some of the 'selected layers' are not in the 'layers' vector.")
   
   g.list <- lapply(sel.layers,
-                   function(x) induced_subgraph(graph, V(graph)[layers == x] ))
+                   function(x) igraph::induced_subgraph(graph, V(graph)[layers == x] ))
   names(g.list) <- sel.layers
   g.list
 }
diff --git a/R/frag.get.layers.pair.R b/R/frag.get.layers.pair.R
index 63a2e10..ea5c6a4 100644
--- a/R/frag.get.layers.pair.R
+++ b/R/frag.get.layers.pair.R
@@ -1,44 +1,44 @@
 
-#'	@export
+
 frag.get.layers.pair  <- function(graph, layer.attr, sel.layers, size.mini=2, mixed.components.only=FALSE)
 {
-  if( ! is.igraph(graph))  stop("Not a graph object")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
+  
   if(! is.logical(mixed.components.only)) stop("The 'mixed.components.only' parameter requires a logical value.")
   if(! is.numeric(size.mini)) stop("The 'size.mini' parameter requires a numerical value.")
   
-  V(graph)$tmp <- vertex_attr(graph, layer.attr)
+  igraph::V(graph)$tmp <- igraph::vertex_attr(graph, layer.attr)
   
-  if(sum(sel.layers %in% V(graph)$tmp) != 2 ){
+  if(sum(sel.layers %in% igraph::V(graph)$tmp) != 2 ){
     stop(paste(c("The values '", sel.layers[1], "' and/or '", sel.layers[2], "' are missing in the '",  layer.attr, "' vertices attribute."),
                sep=" ", collapse = ""))
   }
+  # main function:
+  subgraph <- igraph::induced_subgraph(graph,
+                                       igraph::V(graph)[ igraph::V(graph)$tmp %in% sel.layers ])
   
-  subgraph <- induced_subgraph(graph, V(graph)[ V(graph)$tmp %in% sel.layers ])
-  
-  V(subgraph)$membership <- clusters(subgraph)$membership
+  igraph::V(subgraph)$membership <- igraph::components(subgraph)$membership
   
-  g.list <- decompose(subgraph)   
+  g.list <- igraph::decompose(subgraph)   
   
   if(mixed.components.only == TRUE ){
     sel.components <- sapply(g.list, function(x){
-      (length(unique(V(x)$tmp)) != 1)  & (gorder(x) >= size.mini)
+      (length(unique(V(x)$tmp)) != 1)  & (igraph::gorder(x) >= size.mini)
     })
   }else{
-    sel.components <- sapply(g.list, function(x) gorder(x) >= size.mini)
+    sel.components <- sapply(g.list, function(x) igraph::gorder(x) >= size.mini)
   }
   
-  g <- induced_subgraph(subgraph, V(subgraph)[V(subgraph)$membership %in% which(sel.components)])
-  if(gorder(g) == 0){
-    warnings("There is no mixed component between these layers.")
+  g <- igraph::induced_subgraph(subgraph, 
+             igraph::V(subgraph)[igraph::V(subgraph)$membership %in% which(sel.components)])
+  if(igraph::gorder(g) == 0){
+    warning("There is no mixed component between these layers.")
     return(NULL)
   }
   # remove vertex attribute and return result:
-  delete_vertex_attr(g, "tmp")
+  igraph::delete_vertex_attr(g, "tmp")
 }
 
 
diff --git a/R/frag.get.parameters.R b/R/frag.get.parameters.R
index 2321090..c999525 100644
--- a/R/frag.get.parameters.R
+++ b/R/frag.get.parameters.R
@@ -1,64 +1,67 @@
 
 
 frag.get.parameters <- function(graph, layer.attr){
-  # initial tests:
-  if(! is.igraph(graph)) stop("Not a graph object")
-  if(is_directed(graph)) stop("The 'graph' parameter requires an undirected igraph object.")
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
   
   # retrieve and format 'the layer' attribute:
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
-  V(graph)$layer <- vertex_attr(graph, layer.attr)
+  igraph::V(graph)$layer <- igraph::vertex_attr(graph, layer.attr)
   
   # add edge weight attribute is absent (to avoid an issue with the as_graphnel function):
-  if(is.null(edge_attr(graph, "weight"))){
-    E(graph)$weight <- 1
+  if(is.null(igraph::edge_attr(graph, "weight"))){
+    igraph::E(graph)$weight <- 1
   }
   
   # test of there are two layers:
-  if(length(unique(V(graph)$layer)) != 2) warning("The graph does not have two layers, disturbance and balance values will be meaningless.")
+  if(length(unique(igraph::V(graph)$layer)) != 2) warning("The graph does not have two layers, disturbance and balance values will be meaningless.")
   
   # balance: proportion of non-disturbed pieces in the two layers:
-  v1 <- V(graph)[V(graph)$layer == unique(V(graph)$layer)[1]]
-  v2 <- V(graph)[V(graph)$layer == unique(V(graph)$layer)[2]]
-  subgraph <- subgraph.edges(graph, E(graph)[ ! v1 %--% v2 ])
-  balance <- (table(V(subgraph)$layer) / sum(table(V(subgraph)$layer)) )[1]
+  v1 <- igraph::V(graph)[igraph::V(graph)$layer == unique(igraph::V(graph)$layer)[1]]
+  v2 <- igraph::V(graph)[igraph::V(graph)$layer == unique(igraph::V(graph)$layer)[2]]
+  subgraph <- igraph::subgraph.edges(graph, igraph::E(graph)[ ! v1 %--% v2 ])
+  balance <- (table(igraph::V(subgraph)$layer) / sum(table(igraph::V(subgraph)$layer)) )[1]
   balance <- round(balance, 2)
   
   # components balance:
-  compo.balance <- sapply(decompose(subgraph), function(x) V(x)$layer[1])
+  compo.balance <- sapply(igraph::decompose(subgraph), 
+                          function(x) igraph::V(x)$layer[1])
   compo.balance <- round(table(compo.balance)[1] / sum(table(compo.balance)), 2) 
   
   # disturbance: number of pieces which might have move:
   g.list <- frag.get.layers.pair(graph, "layer", unique(V(graph)$layer), mixed.components.only = TRUE)
   disturbance <- 0
   if(! is.null(g.list)){
-    g.list <- decompose(g.list)
+    g.list <- igraph::decompose(g.list)
     g.list <- sapply(g.list, function(x)
-      table(factor(V(x)$layer, levels = unique(V(graph)$layer))) )
+      table(factor(igraph::V(x)$layer, levels = unique(igraph::V(graph)$layer))) )
     # replace the count of the more represented layer in each component by NA:
     g.list <- apply(g.list, 2, function(x){ x[order(x)][2] <- NA ; x })
     # sum of the count for the less represented layer in each component:
-    disturbance <- sum(g.list, na.rm = TRUE) / gorder(graph)
+    disturbance <- sum(g.list, na.rm = TRUE) / igraph::gorder(graph)
     disturbance <- round(disturbance, 2)
   } 
   
   # degree of aggregation of the edges on the components:
-  aggreg.factor <- 1 - 1/(1 + sd(sapply(decompose(graph), gsize)))
+  aggreg.factor <- 1 - 1/(1 + stats::sd(sapply(igraph::decompose(graph), igraph::gsize)))
   aggreg.factor <- round(aggreg.factor, 2)
   
-  res <- list(n.components = clusters(graph)$no,
-              vertices =  gorder(graph),
-              edges = gsize(graph),
-              balance = balance,
-              components.balance = compo.balance,
-              disturbance = disturbance, 
-              aggreg.factor = aggreg.factor,
-              planar =  boyerMyrvoldPlanarityTest(as_graphnel(graph)))
-  
+  # planarity (if the RBGL package is installed)
+  if (requireNamespace("RBGL", quietly = TRUE)) {
+    is.planar <- RBGL::boyerMyrvoldPlanarityTest(igraph::as_graphnel(graph))
+  }else{
+    warning("The RBGL package is not installed, the `planarity` value is indeterminated and set to FALSE.")
+    is.planar <- FALSE
+  }
+  # list results:
+  res <- list("n.components" = igraph::components(graph)$no,
+              "vertices" =  igraph::gorder(graph),
+              "edges" = igraph::gsize(graph),
+              "balance" = balance,
+              "components.balance" = compo.balance,
+              "disturbance" = disturbance, 
+              "aggreg.factor" = aggreg.factor,
+              "planar" = is.planar)
   # format and return results:
   lapply(res, c, use.names = FALSE)
 }
diff --git a/R/frag.graph.plot.R b/R/frag.graph.plot.R
index d3d045f..7301bc5 100644
--- a/R/frag.graph.plot.R
+++ b/R/frag.graph.plot.R
@@ -1,42 +1,41 @@
 frag.graph.plot  <- function(graph, layer.attr,  ...){
-    if(! is.igraph(graph) )  stop("Not a graph object")
-    if(is.null(graph_attr(graph, "frag_type"))) stop("The 'frag_type' graph attribute is missing")
-    if(is.null(layer.attr)) stop("'layer.attr' is required")
-    if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-    if(is.null(vertex_attr(graph, layer.attr))){
-      stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-    } 
-  
-    V(graph)$layers <- vertex_attr(graph, layer.attr)
+   # tests:
+   .check.frag.graph(graph)
+   .check.layer.argument(graph, layer.attr)
+    if(is.null(igraph::graph_attr(graph, "frag_type"))) stop("The 'frag_type' graph attribute is missing")
+   
+   # main function:
+    igraph::V(graph)$layers <- igraph::vertex_attr(graph, layer.attr)
     
-    nLayers <- length(unique(V(graph)$layers))
-    colors <- c("firebrick2","darkorchid4", "chartreuse3", "darkorange3", "brown3",
-                "darkgoldenrod2","darkolivegreen3", "darksalmon", rainbow(nLayers))
+    nLayers <- length(unique(igraph::V(graph)$layers))
+    colors <- c("firebrick2", "darkorchid4", "chartreuse3", "darkorange3", "brown3",
+                "darkgoldenrod2","darkolivegreen3", "darksalmon")
     # default edge color:
-    E(graph)$color <- "grey"
+    igraph::E(graph)$color <- "grey"
     
-    layers <- sort(unique(V(graph)$layers))
+    layers <- sort(unique(igraph::V(graph)$layers))
     
-    if(graph_attr(graph, "frag_type") == "connection and similarity relations"){
-        graph <- add_layout_(graph, with_fr(), component_wise())
-        E(graph)$color <- as.character(factor(E(graph)$type_relation, labels = c("green", "gray")))
-    } else if(graph_attr(graph, "frag_type") == "similarity relations"){
-      graph <- add_layout_(graph, with_fr(), component_wise())
+    if(igraph::graph_attr(graph, "frag_type") == "connection and similarity relations"){
+        graph <- igraph::add_layout_(graph, igraph::with_fr(), igraph::component_wise())
+        igraph::E(graph)$color <- as.character(factor(igraph::E(graph)$type_relation, labels = c("green", "gray")))
+    } else if(igraph::graph_attr(graph, "frag_type") == "similarity relations"){
+      graph <- igraph::add_layout_(graph, igraph::with_fr(), igraph::component_wise())
     } else if(length(layers) == 2){ 
       # prepare coordinates if the graph has two layers:
-      coords <- data.frame(layer = V(graph)$layers, miny = 0, maxy = 100) 
+      coords <- data.frame(layer = igraph::V(graph)$layers, miny = 0, maxy = 100) 
       coords[coords$layer == layers[1],]$miny <- 51 
       coords[coords$layer == layers[2],]$maxy <- 49 
-      graph$layout <- layout_with_fr(graph,  niter= 1000, weights =  NULL,
+      graph$layout <- igraph::layout_with_fr(graph,  niter= 1000, weights =  NULL,
                                      miny = coords$miny, maxy = coords$maxy  )
     }
     
-    plot(graph, 
-         vertex.color = as.character(factor(V(graph)$layers, labels = colors[1:nLayers] )),  
+    igraph::plot.igraph(graph, 
+         vertex.color = as.character(factor(igraph::V(graph)$layers,
+                                            labels = colors[1:nLayers] )),  
          vertex.label = NA, 
          vertex.size = 4.5,
          edge.width = 2,
-         edge.color = E(graph)$color,
+         edge.color = igraph::E(graph)$color,
          ...)
     invisible(NULL)
 }
diff --git a/R/frag.layers.admixture.R b/R/frag.layers.admixture.R
index ebeae28..286587a 100644
--- a/R/frag.layers.admixture.R
+++ b/R/frag.layers.admixture.R
@@ -1,26 +1,23 @@
 frag.layers.admixture <- function(graph, layer.attr){
   # output : value [0;1]. 0 = "unmixed layers", 1 = "highly mixed layers"
-  if(! is.igraph(graph))  stop("Not a graph object")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
   
   # extract the user-defined layer attribute and reintegrate it as a vertices attribute named "layer":
-  layers <- vertex_attr(graph, layer.attr)
-  V(graph)$layer <- layers
+  layers <- igraph::vertex_attr(graph, layer.attr)
+  igraph::V(graph)$layer <- layers
   layers <- unique(layers)
   
   # Conditionnal tests in function of the number of layers:
   if(length(layers) < 2) stop("At least two different layers are required.")
   
   if(length(layers) == 2){
-    if(is.null(E(graph)$weight)) stop("The edges must be weighted (using the 'frag.edges.weighting' function).")
+    if(is.null(igraph::E(graph)$weight)) stop("The edges must be weighted (using the 'frag.edges.weighting' function).")
     results <- c(admixture = 1 - sum(frag.layers.cohesion(graph, "layer")))
     return(results)
   } else { # if length(layers) > 2
-    pairs <- combn(layers, 2) 
+    pairs <- utils::combn(layers, 2) 
     warning("More than 2 layers: the 'frag.edges.weighting' function has been applied to each pair of layers.")
     results <- sapply(1:ncol(pairs), function(x){
       gsub <- frag.get.layers.pair(graph, layer.attr, c(pairs[1, x], pairs[2, x]))
diff --git a/R/frag.layers.cohesion.R b/R/frag.layers.cohesion.R
index 4eabd83..edba6cc 100644
--- a/R/frag.layers.cohesion.R
+++ b/R/frag.layers.cohesion.R
@@ -1,45 +1,42 @@
 .cohesion.for.two.layers <- function(g, layers){
-  v1 <- V(g)$layers == layers[1]
-  v2 <- V(g)$layers == layers[2]
+  v1 <- igraph::V(g)$layers == layers[1]
+  v2 <- igraph::V(g)$layers == layers[2]
   
   # extract three subgraphs:
-  g1 <-  subgraph.edges(g, E(g)[ V(g)[v1] %--% V(g)[v1] ])
-  g2 <-  subgraph.edges(g, E(g)[ V(g)[v2] %--% V(g)[v2] ])
-  g12 <- subgraph.edges(g, E(g)[ V(g)[v1] %--% V(g)[v2] ])
+  g1 <-  igraph::subgraph.edges(g, igraph::E(g)[ igraph::V(g)[v1] %--% igraph::V(g)[v1] ])
+  g2 <-  igraph::subgraph.edges(g, igraph::E(g)[ igraph::V(g)[v2] %--% igraph::V(g)[v2] ])
+  g12 <- igraph::subgraph.edges(g, igraph::E(g)[ igraph::V(g)[v1] %--% igraph::V(g)[v2] ])
   
   # compute the cohesion values:
-  res1 <-  sum(v1, E(g1)$weight) / sum(gorder(g), E(g)$weight)
-  res2 <-  sum(v2, E(g2)$weight) / sum(gorder(g), E(g)$weight)
+  res1 <-  sum(v1, igraph::E(g1)$weight) / sum(igraph::gorder(g), igraph::E(g)$weight)
+  res2 <-  sum(v2, igraph::E(g2)$weight) / sum(igraph::gorder(g), igraph::E(g)$weight)
   
-  # format and return the results:
+  # format and return results:
   c(res1, res2)
 }
 
 frag.layers.cohesion <- function(graph, layer.attr){
   # output : value [0;1].
-  if(! is.igraph(graph)) stop("Not a graph object")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute.", sep=""))
-  }
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
   
   # delete singletons:
-  graph <- delete_vertices(graph, degree(graph) == 0)
+  graph <- igraph::delete_vertices(graph, degree(graph) == 0)
   
   # extract the user-defined layer attribute and reintegrate it as a vertices attribute named "layer":
-  layers <- vertex_attr(graph, name = layer.attr)
-  V(graph)$layers <- layers
+  layers <- igraph::vertex_attr(graph, name = layer.attr)
+  igraph::V(graph)$layers <- layers
   layers <- sort(unique(layers))
   # Conditionnal tests in function of the number of layers:
   if(length(layers) < 2){
     warning("At least two different layers are required.")
     return(c(NA, NA))
   }
-  pairs <- combn(layers, 2) 
+  pairs <- utils::combn(layers, 2) 
   
   if(length(layers) == 2){
-    if(is.null(E(graph)$weight)) stop("The edges must be weighted (using the 'frag.edges.weighting' function).")
+    if(is.null(igraph::E(graph)$weight)) stop("The edges must be weighted (using the 'frag.edges.weighting' function).")
     results <- .cohesion.for.two.layers(graph, layers)
     results <- matrix(results)
   } else{ # if length(layers) > 2
@@ -47,7 +44,7 @@ frag.layers.cohesion <- function(graph, layer.attr){
     results <- sapply(1:ncol(pairs), function(x){
       gsub <- frag.get.layers.pair(graph, layer.attr, c(pairs[1, x], pairs[2, x]))
       gsub <- frag.edges.weighting(gsub, layer.attr)
-      .cohesion.for.two.layers(gsub, unique(V(gsub)$layers))
+      .cohesion.for.two.layers(gsub, unique(igraph::V(gsub)$layers))
     })
   }
   rownames(results) <- c("cohesion1", "cohesion2")
diff --git a/R/frag.observer.failure.R b/R/frag.observer.failure.R
index 3e54e6c..ea03f06 100644
--- a/R/frag.observer.failure.R
+++ b/R/frag.observer.failure.R
@@ -1,12 +1,14 @@
 .randomly.delete.edges <- function(graph, value){ 
   set.seed(10)  # reset the random value to ensure that the same series of edges is remove
-  n.edges <- round(gsize(graph) * value)
-  delete_edges(graph, sample(E(graph), n.edges))
+  n.edges <- round(igraph::gsize(graph) * value)
+  igraph::delete_edges(graph, sample(igraph::E(graph), n.edges))
 }
 
 frag.observer.failure <- function(graph, likelihood){
   # output: a liste of altered graphs with edges removed
-  if(! is.igraph(graph))  stop("Not a graph object")
+  # tests:
+  .check.frag.graph(graph)
+  # main function:  
   if(! is.numeric(likelihood)){
     stop("The 'likelood' parameter requires a numerical value.")
   } else  if(sum(likelihood < 0, likelihood > 1)){
diff --git a/R/frag.path.lengths.R b/R/frag.path.lengths.R
index c75aa40..5dd1195 100644
--- a/R/frag.path.lengths.R
+++ b/R/frag.path.lengths.R
@@ -1,8 +1,10 @@
 frag.path.lengths <- function(graph, cumulative=FALSE){
-  if(! is.igraph(graph)) stop("Not a graph object")
+  # tests:
+  .check.frag.graph(graph)
+  # main function:  
   if(! is.logical(cumulative)) stop("The 'cumulative' parameter requires a logical value.")
   
-  path.vector <- distance_table(graph, directed = FALSE)$res
+  path.vector <- igraph::distance_table(graph, directed = FALSE)$res
   
   if(cumulative){
     path.vector <- path.vector / max(path.vector)
diff --git a/R/frag.relations.by.layers.R b/R/frag.relations.by.layers.R
index ac3d038..e2c7c27 100644
--- a/R/frag.relations.by.layers.R
+++ b/R/frag.relations.by.layers.R
@@ -1,18 +1,17 @@
 frag.relations.by.layers <- function(graph, layer.attr){
-  if(! is.igraph(graph)) stop("Not an igraph object")
-  if(is.null(vertex_attr(graph, layer.attr))) stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute", sep=""))
-  }
-  layers <- vertex_attr(graph, layer.attr)
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
+  
+  # retrieve the layer argument:
+  layers <- igraph::vertex_attr(graph, layer.attr)
   
-  if(is.null(V(graph)$name)){
-    V(graph)$name <- 1:gorder(graph)
+  if(is.null(igraph::V(graph)$name)){
+    igraph::V(graph)$name <- 1:igraph::gorder(graph)
   }
   
-  v.list <- data.frame(v = V(graph)$name, layer = layers)
-  e.list <- data.frame(as_edgelist(graph))
+  v.list <- data.frame("v" = igraph::V(graph)$name, "layer" = layers)
+  e.list <- data.frame(igraph::as_edgelist(graph))
   e.list <- merge(e.list, v.list, by.x = "X2", by.y = "v")
   e.list <- merge(e.list, v.list, by.x = "X1", by.y = "v")
   
diff --git a/R/frag.simul.compare.R b/R/frag.simul.compare.R
index 0759a71..608987d 100644
--- a/R/frag.simul.compare.R
+++ b/R/frag.simul.compare.R
@@ -1,17 +1,17 @@
 
-.run.simul <- function(iter, observed.graph, layer.attr, initial.layer){
+.run.simul <- function(iter, observed.graph, layer.attr, initial.layer, ...){
   res <- lapply(1:iter, function(i){
     g <- frag.simul.process(initial.layers = initial.layer,
                             from.observed.graph = observed.graph,
-                            observed.layer.attr = layer.attr)
+                            observed.layer.attr = layer.attr, ...)
     g <- frag.edges.weighting(g, layer.attr)
     # measure the properties of the graph:
-    inter.layer.e <- E(g)[V(g)[V(g)$layer == 1] %--% V(g)[V(g)$layer == 2]]
+    inter.layer.e <- igraph::E(g)[ igraph::V(g)[igraph::V(g)$layer == 1] %--% igraph::V(g)[igraph::V(g)$layer == 2]]
     c(
-      "edges" = gsize(g),
-      "weightsum" = sum(E(g)$weight),
-      "balance" = c(sort(table(V(g)$layer))[1] / sum(table(V(g)$layer)), use.names=FALSE),
-      "disturbance" = length(inter.layer.e) / gsize(g),
+      "edges" = igraph::gsize(g),
+      "weightsum" = sum(igraph::E(g)$weight),
+      "balance" = c(sort(table(igraph::V(g)$layer))[1] / sum(table(igraph::V(g)$layer)), use.names=FALSE),
+      "disturbance" = length(inter.layer.e) / igraph::gsize(g),
       frag.layers.admixture(g, "layer"),
       "cohesion" = rbind(frag.layers.cohesion(g, "layer"))
     )
@@ -20,27 +20,25 @@
   data.frame(res)
 }  
 
-frag.simul.compare <- function(graph, layer.attr, iter, summarise=TRUE){
-  if(! is.igraph(graph)) stop("Not a graph object")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute.", sep=""))
-  }
+frag.simul.compare <- function(graph, layer.attr, iter, summarise=TRUE, ...){
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
+  
   if(iter < 30) stop("At least 30 iterations are required.")
   if(! is.logical(summarise)){
     stop("A logical value is required for the 'summary' parameter.")
   }
-  
-  resH1 <- .run.simul(iter, graph, layer.attr, 1)
-  resH2 <- .run.simul(iter, graph, layer.attr, 2)
+  # main funtion:
+  resH1 <- .run.simul(iter, graph, layer.attr, 1, ...)
+  resH2 <- .run.simul(iter, graph, layer.attr, 2, ...)
   
   if(! summarise){
-    return(list(h1.data = resH1, h2.data = resH2))
+    return(list("h1.data" = resH1, "h2.data" = resH2))
   }
   stats <- frag.simul.summarise(graph, layer.attr, resH1, resH2)
   # print the summary table:
   print(stats)
   #return silently the results as a list:
-  invisible(list(h1.data = resH1, h2.data = resH2, summary = stats))
+  invisible(list("h1.data" = resH1, "h2.data" = resH2, "summary" = stats))
 }
diff --git a/R/frag.simul.process.R b/R/frag.simul.process.R
index 7b56037..5ef8aab 100644
--- a/R/frag.simul.process.R
+++ b/R/frag.simul.process.R
@@ -1,7 +1,7 @@
 .subsetsum <- function(x, target, i = 1){
   while(i != length(x)){
     s <- sum(x[1:i], na.rm = T)
-    if(s==target) break
+    if(s == target) break
     if(s > target) x[i] <- NA
     i <- i + 1
   }
@@ -11,34 +11,30 @@
 
 .select.component <- function(g, aggreg.factor){
   # high aggreg.factor favors bigger components:
-  proba <- 1/ (1 + (1:clusters(g)$no * aggreg.factor))
-  sample(order(clusters(g)$csize, decreasing = T), 1, prob = proba)
-}
-
-.connect.neighbors.if.planar <- function(g, v, v.to.add.name, neighbors){
-  g.tmp <- induced_subgraph(g, V(g)[ V(g)$object.id == V(g)[v]$object.id ])
-  g.tmp <- add_vertices(g.tmp, 1, attr = list(name = v.to.add.name))
-  g.tmp <- add_edges(g.tmp, as.character(c(rbind(v.to.add.name, neighbors)) ))
-  
-  if(boyerMyrvoldPlanarityTest(as_graphnel(g.tmp))){
-    g <- add_edges(g, c(rbind(v.to.add.name, neighbors)) )
-  }
-  g    
+  proba <- 1/ (1 + (1:igraph::components(g)$no * aggreg.factor))
+  sample(order(igraph::components(g)$csize, decreasing = T), 1, prob = proba)
 }
 
 .connect.neighbors <- function(g, v, v.to.add.name, edges, planar){
+  # if there are other neighbours than the new vertex, try to connect one of them with the new vertex.
   # identify the neighbors of the target vertex:
-  neighbors <- neighborhood(g, order = 1, nodes = v, mindist=1)
+  neighbors <- igraph::ego(g, order = 1, nodes = v, mindist = 1)
   neighbors <- unlist(neighbors)
   neighbors <- neighbors[ ! neighbors == v.to.add.name]
   neighbors <- as.character(neighbors)
-  neighbors <- sample(neighbors, sample(0:length(neighbors), 1) )
-  # if there are other neighbours than the new vertex, try connection:
-  if(length(neighbors) > 0 & (gsize(g) + length(neighbors)) < edges ){
+  # select none or one of the neighbors:
+  neighbors <- sample(neighbors, sample(0:length(neighbors), 1))
+  if(length(neighbors) > 0 & (igraph::gsize(g) + 1) < edges ){
     if(planar){
-      g <- .connect.neighbors.if.planar(g, v, v.to.add.name, neighbors)
-    }else{
-      g <- add_edges(g, c(rbind(v.to.add.name, neighbors)) ) # planarity test
+      # make a temporary version of the modified component of the graph:
+      g.tmp <- igraph::induced_subgraph(g, igraph::V(g)[ igraph::V(g)$object.id == igraph::V(g)[v]$object.id ])
+      g.tmp <- igraph::add_edges(g.tmp, as.character(c(rbind(v.to.add.name, neighbors)) ))
+      # if the modified component is still planar, then connect the new vertex and the neighbor: 
+      if(RBGL::boyerMyrvoldPlanarityTest(igraph::as_graphnel(g.tmp))){
+        g <- igraph::add_edges(g, as.character(c(rbind(v.to.add.name, neighbors)) ))
+      }
+    }else{ # if no planarity constraint:
+      g <- igraph::add_edges(g, c(rbind(v.to.add.name, neighbors)) )
     }
   }
   g
@@ -47,14 +43,15 @@
 .add.fragment <- function(g, n.components, edges, connect.neighbors, planar, aggreg.factor){
   # select component:
   compo.to.fragment <- .select.component(g, aggreg.factor)
-  v.to.fragment <- V(g)[V(g)$object.id == compo.to.fragment]$name
+  
+  v.to.fragment <- igraph::V(g)[ igraph::V(g)$object.id == compo.to.fragment]$name
   # select the target vertex in the component:
   v <-  sample(as.character(v.to.fragment), 1)
   # add a new vertex and an edge with the target vertex:
-  v.to.add.name <- gorder(g) + 1
-  g <- add_vertices(g, 1, attr = list(name = v.to.add.name, 
-                                      object.id = V(g)[v]$object.id ))
-  g <- add_edges(g, c(v, v.to.add.name))
+  v.to.add.name <- igraph::gorder(g) + 1
+  g <- igraph::add_vertices(g, 1, attr = list("name" = v.to.add.name, 
+                                      "object.id" = igraph::V(g)[v]$object.id ))
+  g <- igraph::add_edges(g, c(v, v.to.add.name))
   
   # try to connect the new vertex and the neighbors of the target vertex:
   if(connect.neighbors){
@@ -64,23 +61,24 @@
 }
 
 
+
+
 .main <- function(n.components, vertices, edges, balance, disturbance, aggreg.factor, planar){
   # Initialize graph:  
-  g <- n.components * make_graph(c(1, 2), directed = FALSE)
-  V(g)$name <- 1:gorder(g)
-  V(g)$object.id <- clusters(g)$membership
-  
+  g <- n.components * igraph::make_graph(c(1, 2), directed = FALSE)
+  igraph::V(g)$name <- 1:igraph::gorder(g)
+  igraph::V(g)$object.id <- igraph::components(g)$membership
   # Build graph:
   if(is.infinite(vertices)){ # only edge count constraint
     if(vertices/2 < n.components){
       stop("Increase 'vertices' or decrease 'n.components'.")}
-    while(gsize(g) < edges){
+    while(igraph::gsize(g) < edges){
       g <- .add.fragment(g, n.components, edges, connect.neighbors=T, planar, aggreg.factor)
     }
   } else if(is.infinite(edges)){ # only vertices count constraint
     if(edges < n.components){
       stop("Increase 'edges' or decrease 'n.components'.")}
-    while(gorder(g) < vertices){
+    while(igraph::gorder(g) < vertices){
       g <- .add.fragment(g, n.components, edges, connect.neighbors=T, planar, aggreg.factor)
     }
   } else if(! is.infinite(vertices) & ! is.infinite(edges) ){ # vertices & edges constraints
@@ -95,135 +93,138 @@
       stop("Irrelevant parameters, decrease 'edges' or increase 'vertices'.")
     }  
     
-    while(gorder(g) < vertices & gsize(g) < edges ){
+    while(igraph::gorder(g) < vertices & igraph::gsize(g) < edges ){
       g <- .add.fragment(g, n.components, edges, connect.neighbors=F, planar, aggreg.factor)
     }
-    # check if the number of supplementary edges can conserve planarity:
+    # check if the number of supplementary edges conserve planarity:
     e.existing <- 0 # initialize the value
     e.max <- edges  # initialize the value
-    clus <- clusters(g)$csize
-    gsub <- induced_subgraph(g,
-                             V(g)[ clusters(g)$membership %in% which(clus > 2)])
-    gsub <- decompose(gsub)
+    components.size <- igraph::components(g)$csize
+    gsub <- igraph::induced_subgraph(g,
+                             igraph::V(g)[ igraph::components(g)$membership %in% which(components.size > 2)])
+    gsub <- igraph::decompose(gsub)
     if(length(gsub) > 0){
-      e.existing <- sapply(gsub, gsize)
-      e.max <- sapply(clus[clus > 2], function(x)  3*x-6) # e max for planar graphs
+      e.existing <- sapply(gsub, igraph::gsize)
+      e.max <- sapply(components.size[components.size > 2], function(x)  3*x-6) # e max for planar graphs
       e.max <- sum(e.max - e.existing)
     }
     
-    if(edges - gsize(g) > e.max){
-      # message("Few or no solution possible for these parameters. Consider changing them. The graph generated has less edges than the 'edges' parameter.")
-      # edges <- gsize(g) + e.max # the nr of edges is reduced
+    if(edges - igraph::gsize(g) > e.max){
       stop("No solution for these parameters, given the current random result. Consider removing the edge or vertex constraint.")
     }
     while(gsize(g) < edges){
       # select a component:
       selected.component <- .select.component(g, aggreg.factor)
-      v.to.connect <- V(g)[V(g)$object.id == selected.component]$name
+      v.to.connect <- igraph::V(g)[ igraph::V(g)$object.id == selected.component]$name
       # select two vertices in the component:
       v.to.connect <-  sample(as.character(v.to.connect), 2)
-      if(length(E(g)[v.to.connect[1] %--% v.to.connect[2]]) == 0){# if edge not exists
-        # check planarity before adding edge:
+      if(length(igraph::E(g)[v.to.connect[1] %--% v.to.connect[2]]) == 0){# if the edge does not exist yet
+        # check planarity of the component before adding edge:
         if(planar){
-          g.tmp <- induced_subgraph(g, V(g)[ V(g)$object.id == selected.component])
-          g.tmp <- add_edges(g.tmp, v.to.connect)
+          g.tmp <- igraph::induced_subgraph(g, igraph::V(g)[ igraph::V(g)$object.id == selected.component])
+          g.tmp <- igraph::add_edges(g, v.to.connect)
           
-          if(boyerMyrvoldPlanarityTest(as_graphnel(g.tmp))){
-            g <- add_edges(g, v.to.connect)
+          if(RBGL::boyerMyrvoldPlanarityTest(igraph::as_graphnel(g.tmp))){
+            g <- igraph::add_edges(g, v.to.connect)
           }
         }else{
-          g <- add_edges(g, v.to.connect)
+          g <- igraph::add_edges(g, v.to.connect)
         }
       }
-    }
-    
-  }
+    } # end of the while loop
+  } # end of the elseif "vertices + edges constraints"
   g
 }
 
 
 .add.disturbance <- function(g, nr.v.to.disturb, asymmetric.transport.from){
   # default behaviour:
-  v.to.disturb <- sample(seq(1, gorder(g)), nr.v.to.disturb)
+  v.to.disturb <- sample(seq(1, igraph::gorder(g)), nr.v.to.disturb)
   # if asymmetric.transport.from is set:
   if( ! is.null(asymmetric.transport.from)){
-    if(nr.v.to.disturb <= length(V(g)[V(g)$layer == asymmetric.transport.from]) ){
-      v.to.disturb <- sample(V(g)[V(g)$layer == asymmetric.transport.from], nr.v.to.disturb)
+    if(nr.v.to.disturb <= length(igraph::V(g)[ igraph::V(g)$layer == asymmetric.transport.from]) ){
+      v.to.disturb <- sample(igraph::V(g)[ igraph::V(g)$layer == asymmetric.transport.from], nr.v.to.disturb)
     } else{
       stop("The number of fragments for asymmetric transport exceeds the number of fragments in this layer. Reduce the 'disturbance' value.")
     }
   } 
   # reverse layer values of the selected vertices:
-  V(g)[v.to.disturb]$layer <- as.character(factor(V(g)[v.to.disturb]$layer,
+  igraph::V(g)[v.to.disturb]$layer <- as.character(factor(V(g)[v.to.disturb]$layer,
                                                   levels = c(1,2), labels = c(2,1)))
   g
 }
 
-frag.simul.process <- function(initial.layers=2, n.components, vertices=Inf, edges=Inf, balance=.5, components.balance=.5, disturbance=0, aggreg.factor=0, planar=TRUE, asymmetric.transport.from=NULL, from.observed.graph=NULL, observed.layer.attr=NULL){
+
+ 
+
+frag.simul.process <- function(initial.layers=2, n.components, vertices=Inf, edges=Inf, balance=.5, components.balance=.5, disturbance=0, aggreg.factor=0, planar=FALSE, asymmetric.transport.from=NULL, from.observed.graph=NULL, observed.layer.attr=NULL){
   
-  #  If requested input parameters from observed graph (except the number of edges):
+  if(! is.logical(planar)) stop("The 'planar' argument must be logical.")
+  if(planar==TRUE & (! requireNamespace("RBGL", quietly=TRUE))){
+    stop("To use the `planar` constraint, the RBGL package is required.")
+  }
+  
+  #  If required by the user, use parameters from the observed graph (except the number of edges):
   if( ! is.null(from.observed.graph) & ! is.null(observed.layer.attr)){
     if( ! is.character(observed.layer.attr))  stop("The parameter 'observed.layer.attr' requires a character value.")
-    if( ! observed.layer.attr %in% names(vertex_attr(from.observed.graph)) ){
-      stop(paste("No '", observed.layer.attr, "' vertices attribute", sep=""))
+    if( ! observed.layer.attr %in% names(igraph::vertex_attr(from.observed.graph)) ){
+      stop(paste("No '", observed.layer.attr, "' vertices attribute.", sep=""))
     }
     if(length(unique(vertex_attr(from.observed.graph, observed.layer.attr))) != 2){
-      stop("The layer attribute of the observed graph must contain two layers.")
+      stop("The `layer` attribute of the observed graph must contain two layers.")
     }
-    # run the get.parameters function:
+    # retrieve the observed graph's values:
     params <- frag.get.parameters(from.observed.graph, observed.layer.attr)
-    # input the observed parameters:
-    n.components <- params$n.components
-    vertices <- params$vertices
-    balance <- params$balance
-    components.balance <- params$components.balance
-    disturbance <- params$disturbance
-    aggreg.factor <- params$aggreg.factor
-    planar <- params$planar
+    # set the parameters with observed values if not already set by the user:
+    if(missing(n.components)) n.components <- params$n.components
+    if(missing(vertices)) vertices <- params$vertices
+    if(missing(balance)) balance <- params$balance
+    if(missing(components.balance)) components.balance <- params$components.balance
+    if(missing(disturbance)) disturbance <- params$disturbance
+    if(missing(aggreg.factor)) aggreg.factor <- params$aggreg.factor
+    if(missing(planar))  planar <- params$planar
   }
   
   # BEGIN Tests:
-  if(! is.logical(planar)) stop("The 'planar' argument must be logical.")
   if(is.null(n.components)) stop("The 'n.components' parameter is required.")
-  
+
   if(! is.numeric(balance)){
     stop("The 'balance' argument requires a numerical value.")
   } else if(balance <= 0 | balance >= 1){
     stop("'balance' values must range in ]0;1[")
   }
-  
+
   if(! is.numeric(components.balance)){
     stop("The 'components.balance' argument requires a numerical value.")
   } else if(components.balance <= 0 | components.balance >= 1){
     stop("'components.balance' values must range in ]0;1[")
   }
-  
+
   if(! is.numeric(disturbance)){
     stop("The 'disturbance' argument requires a numerical value.")
   } else if(disturbance < 0 | disturbance > 1){
     stop("'disturbance' values must range in [0;1].")
   }
-  
+
   if(is.infinite(vertices) & is.infinite(edges)){
     stop("At least one of the parameters 'vertices' or 'edges' is required.")
   }
   if(! initial.layers %in% c(1, 2)){
     stop("The 'initial.layers' parameter requires a numerical value of 1 or 2.")
   }
-  
+
   if(! is.numeric(aggreg.factor)){
     stop("The 'disturbance' argument requires a numerical value.")
   } else if(aggreg.factor > 1 | aggreg.factor < 0 ){
     stop("The 'aggreg.factor' parameter must range in [0;1].")
   }
-  
+
   if( ! is.null(asymmetric.transport.from) ){
     if(! asymmetric.transport.from %in% c(1, 2, "1", "2")){
       stop("The 'asymmetric.transport.from' parameter must have a value in 1 or 2.")
-    } 
+    }
   }
-  
-  # END tests.
+  # END tests
   
   # BEGIN main body of the function:
   
@@ -231,20 +232,20 @@ frag.simul.process <- function(initial.layers=2, n.components, vertices=Inf, edg
     g <- .main(n.components, vertices, edges, balance, disturbance, aggreg.factor, planar) 
     
     # BALANCE. Determine layer size:
-    v.layer1 <- round(gorder(g) * balance)
+    v.layer1 <- round(igraph::gorder(g) * balance)
     
     # search possible combinations of components and use the first one:
-    sel.clusters <- clusters(g)$csize
-    names(sel.clusters) <- seq_len(length(sel.clusters))
-    sel.clusters <- .subsetsum(sample(sel.clusters), v.layer1) # randomize order
-    sel.clusters <- names(sel.clusters)
+    sel.components <- igraph::components(g)$csize
+    names(sel.components) <- seq_len(length(sel.components))
+    sel.components <- .subsetsum(sample(sel.components), v.layer1) # randomize order
+    sel.components <- names(sel.components)
     
     # assign layers:
-    V(g)$layer <- 2
-    V(g)[ V(g)$object.id %in% sel.clusters ]$layer <- 1
+    igraph::V(g)$layer <- 2
+    igraph::V(g)[ igraph::V(g)$object.id %in% sel.components ]$layer <- 1
     
     # ADD DISTURBANCE:
-    nr.v.to.disturb <- round(gorder(g) * disturbance)
+    nr.v.to.disturb <- round(igraph::gorder(g) * disturbance)
     if(nr.v.to.disturb > 0){
       g <- .add.disturbance(g, nr.v.to.disturb, asymmetric.transport.from)
     }
@@ -277,18 +278,19 @@ frag.simul.process <- function(initial.layers=2, n.components, vertices=Inf, edg
                       aggreg.factor,
                       planar) 
     # mark and merge the two graphs:
-    V(g.layer1)$layer <- 1
-    V(g.layer2)$layer <- 2
-    V(g.layer2)$name <- paste(V(g.layer2)$name, ".2", sep="")
-    g <- g.layer1 %du% g.layer2
+    igraph::V(g.layer1)$layer <- 1
+    igraph::V(g.layer2)$layer <- 2
+    igraph::V(g.layer2)$name <- paste(igraph::V(g.layer2)$name, ".2", sep="")
+    g <- igraph::disjoint_union(g.layer1, g.layer2)
     # ADD DISTURBANCE:
-    nr.v.to.disturb <- round(gorder(g) * disturbance)
+    nr.v.to.disturb <- round(igraph::gorder(g) * disturbance)
     if(nr.v.to.disturb > 0){
       g <- .add.disturbance(g, nr.v.to.disturb, asymmetric.transport.from)
     }
   }
+  # finalize and return the graph:
   g <- frag.edges.weighting(g, layer.attr="layer")
-  g <- delete_vertex_attr(g, "which")
+  g <- igraph::delete_vertex_attr(g, "which")
   g$frag_type <- "cr"
   g
 }
diff --git a/R/frag.simul.summarise.R b/R/frag.simul.summarise.R
index df468a7..402c7e6 100644
--- a/R/frag.simul.summarise.R
+++ b/R/frag.simul.summarise.R
@@ -4,7 +4,7 @@
   h1.values <- unlist(h1.values)
   h2.values <- unlist(h2.values)
   # difference between H1 and H2
-  wilcox.res <- wilcox.test(h1.values, h2.values, exact=FALSE)$p.value
+  wilcox.res <- stats::wilcox.test(h1.values, h2.values, exact=FALSE)$p.value
   
   if(is.nan(wilcox.res)) return(c(NA,NA,NA,NA))
   
@@ -38,12 +38,10 @@ frag.simul.summarise <- function(graph, layer.attr, res.h1, res.h2,
                                  cohesion2.attr = "cohesion2",
                                  admixture.attr = "admixture"){
   # todo: add params:  
-  if(! is.igraph(graph)) stop("Not a graph object")
-  if(is.null(vertex_attr(graph, layer.attr)))   stop("'layer.attr' is missing or does not correspond to a vertex attribute of the graph.")
-  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
-  if( ! layer.attr %in% names(vertex_attr(graph)) ){
-    stop(paste("No '", layer.attr, "' vertices attribute.", sep=""))
-  }
+  # tests:
+  .check.frag.graph(graph)
+  .check.layer.argument(graph, layer.attr)
+  
   if(is.null(res.h1) | is.null(res.h2)){
     stop("res.h1 and res.h2 are mandatory.")
   }
@@ -70,7 +68,7 @@ frag.simul.summarise <- function(graph, layer.attr, res.h1, res.h2,
   obs.params <- c(frag.get.parameters(graph, layer.attr),
                   frag.layers.admixture(graph, layer.attr),
                   "cohesion" = frag.layers.cohesion(graph, layer.attr),
-                  "weightsum" = sum(E(graph)$weight))
+                  "weightsum" = sum(igraph::E(graph)$weight))
   
   if(sum(! colnames(res.h1) %in% names(obs.params)) != 0){
     warning("Some simulated parameters are missing in the observed graph.")
diff --git a/R/make_frag_object.R b/R/make_frag_object.R
index 52b17fe..97b3094 100644
--- a/R/make_frag_object.R
+++ b/R/make_frag_object.R
@@ -18,24 +18,24 @@ setGeneric(
 
 make_frag_set_validation <- function(object)
 {
-  if( ncol(object@cr.df) == 1 & ncol(object@sr.df) == 1 ){
+  if( ncol(object@df.cr) == 1 & ncol(object@df.sr) == 1 ){
     stop("At least one of the 'cr' and 'sr' arguments is required.")
   }
   
   if(  object@frag_type %in% c("sr", "crsr")) {
-    if( ! ncol(object@sr.df) > 1){
+    if( ! ncol(object@df.sr) > 1){
       stop("The data frame for the 'sr' argument must have at least two columns.")
     }
-    if( sum( ! object@sr.df[,1] %in% object@fragments.df[, 1]) != 0){
+    if( sum( ! object@df.sr[,1] %in% object@fragments.df[, 1]) != 0){
       stop("Some objects in 'sr' are not documented in the 'fragments' data frame.")
     }
   }  
 
   if( object@frag_type %in% c("cr", "crsr") ) {
-    if( ! ncol(object@cr.df) > 1){
+    if( ! ncol(object@df.cr) > 1){
       stop("The data frame for the 'cr' argument must have at least two columns.")
     }
-    if( sum( ! c(object@cr.df[,1], object@cr.df[,2]) %in% object@fragments.df[, 1]) != 0){
+    if( sum( ! c(object@df.cr[,1], object@df.cr[,2]) %in% object@fragments.df[, 1]) != 0){
       stop("Some objects in 'cr' are not documented in the 'fragments' data frame.")
     }
   } 
@@ -49,8 +49,8 @@ make_frag_set_validation <- function(object)
 setClass(
   Class="Frag.object", 
   representation=representation(
-    cr.df="matrix",
-    sr.df="matrix",
+    df.cr="matrix",
+    df.sr="matrix",
     fragments.df="data.frame",
     frag_type="character"
   ),
@@ -66,9 +66,9 @@ setMethod(
       "------ Frag.object ------",
       "\n*   Frag_type = ", object@frag_type,
       "\n*   N fragments = ",
-      length(unique( na.omit(c(object@cr.df[,1],
-                               object@cr.df[,2],
-                               object@sr.df[,1]))) ),
+      length(unique( stats::na.omit(c(object@df.cr[,1],
+                                      object@df.cr[,2],
+                                      object@df.sr[,1]))) ),
       "\n-------------------------",
       sep=""))
   }
@@ -82,10 +82,12 @@ setMethod(
     if( object@frag_type == "sr" ){
       stop("No available data for the connection relationships.")
     }
-    cr.net <- graph_from_data_frame(object@cr.df, directed=FALSE, vertices=object@fragments.df)
-    cr.net <- delete_vertices(cr.net, degree(cr.net, mode="total") == 0)
-    E(cr.net)$type_relation <- "cr"
-    cr.net <- set_graph_attr(cr.net, "frag_type", "connection relations")
+    cr.net <- igraph::graph_from_data_frame(object@df.cr, directed=FALSE,
+                                            vertices=object@fragments.df)
+    cr.net <- igraph::delete_vertices(cr.net, 
+                                      igraph::degree(cr.net, mode="total") == 0)
+    igraph::E(cr.net)$type_relation <- "cr"
+    cr.net <- igraph::set_graph_attr(cr.net, "frag_type", "connection relations")
     return(cr.net)
   }
 )
@@ -99,24 +101,26 @@ setMethod(
       stop("No available data for the similarity relationships.")
     }
     #  'similarity units' ids are recoded to avoid confusion with the fragments ids:
-    object@sr.df[,2] <- as.character(factor(
-      object@sr.df[,2], labels=paste("su", c(1:(0 + length(unique((object@sr.df[,2]))) )))
+    object@df.sr[,2] <- as.character(factor(
+      object@df.sr[,2], labels=paste("su", c(1:(0 + length(unique((object@df.sr[,2]))) )))
     ) )
-    vertices.list <- unique( c(object@sr.df[, 1], object@sr.df[, 2]) ) 
-    sr.net <- simplify(graph_from_data_frame(object@sr.df[, 1:2], directed=TRUE, vertices=vertices.list ))
-    sr.net <- graph_from_adjacency_matrix(bibcoupling(sr.net), diag=FALSE, mode="undirected" )
-    sr.net <- delete_vertices(sr.net, degree(sr.net, mode="total") == 0)
+    vertices.list <- unique( c(object@df.sr[, 1], object@df.sr[, 2]) ) 
+    sr.net <- igraph::simplify(igraph::graph_from_data_frame(object@df.sr[, 1:2], 
+                                                 directed=TRUE, vertices=vertices.list ))
+    sr.net <- igraph::graph_from_adjacency_matrix(igraph::bibcoupling(sr.net),
+                                                  diag=FALSE, mode="undirected" )
+    sr.net <- igraph::delete_vertices(sr.net, igraph::degree(sr.net, mode="total") == 0)
     # vertices attributes:
     fragments.df <- object@fragments.df
     names(fragments.df)[1] <- "name"
     
-    attributes <- merge(                 #retrieve the  graph attributes
-      cbind(name = V(sr.net)$name),
+    attributes <- merge(       #retrieve the  graph attributes
+      cbind("name" = igraph::V(sr.net)$name),
       cbind(fragments.df),
       by="name", sort=FALSE)
-    vertex_attr(sr.net) <- lapply(attributes, as.character) #add vertex attributes
+    igraph::vertex_attr(sr.net) <- lapply(attributes, as.character) #add vertex attributes
     # edge attribute:
-    E(sr.net)$type_relation <- "sr"
+    igraph::E(sr.net)$type_relation <- "sr"
     return(sr.net)
   }
 )
@@ -132,13 +136,14 @@ setMethod(
     cr.net <- make_cr_graph(object)
     sr.net <- make_sr_graph(object)
     
-    crsr.net <- cr.net %u% sr.net
-    crsr.list <- decompose(crsr.net) # get one graph for each component
+    crsr.net <- igraph::union(cr.net, sr.net)
+    crsr.list <- igraph::decompose(crsr.net) # get one graph for each component
     # union of each graph (ie: connection graph) with its complement graph (ie: similarity graph)
     crsr.list <- lapply(crsr.list,
-                            function(x) x %u% complementer(x) )
+                            function(x) igraph::union(x, igraph::complementer(x)) )
     crsr.list <- lapply(crsr.list,
-                            function(x){ set_vertex_attr(x, "name_save", V(x),  V(x)$name )} )
+                            function(x){ igraph::set_vertex_attr(x, "name_save", 
+                                                         igraph::V(x),  igraph::V(x)$name )} )
     # merge all the graphs in the list:
     crsr.net <- Reduce("union", crsr.list)
     
@@ -147,30 +152,30 @@ setMethod(
     fragments.df <- object@fragments.df
     names(fragments.df)[1] <- "name"
     attributes <- merge(                 # retrieve the crsr.net graph attributes
-      cbind(name = V(crsr.net)$name),
+      cbind("name" = igraph::V(crsr.net)$name),
       cbind(fragments.df),
       by="name", sort=FALSE)
-    vertex_attr(crsr.net) <- lapply(attributes, as.character) #add vertex attributes
+    igraph::vertex_attr(crsr.net) <- lapply(attributes, as.character) #add vertex attributes
     
     # 2. edge attributes ####
-    edge.attributes(crsr.net) <- list()  # removing all edges attributes
+    igraph::edge_attr(crsr.net) <- list()  # removing all edges attributes
     # we compare the edges in cr.net and crsr.net and only keep the "connection" edges
     # we got a list of the crsr.net edges also present in cr.net
     rename.edges <- function(edge.list){
       chain <- paste(edge.list[[1]], edge.list[[2]], sep="-")
       return(chain)
     }
-    crsr.net.edgelist <- as_edgelist(crsr.net)
+    crsr.net.edgelist <- igraph::as_edgelist(crsr.net)
     crsr.net.edgelist <- apply(crsr.net.edgelist, 1, rename.edges)
-    cr.net.edgelist <- as_edgelist(cr.net)
+    cr.net.edgelist <- igraph::as_edgelist(cr.net)
     cr.net.edgelist <- apply(cr.net.edgelist, 1, rename.edges)
     
     # setting the "type_relation" edge attribute to "cr"
-    E(crsr.net)$type_relation <- NA
-    E(crsr.net)[ which(crsr.net.edgelist %in% cr.net.edgelist) ]$type_relation <- "cr"
-    E(crsr.net)[ is.na( E(crsr.net)$type_relation ) ]$type_relation <- "sr"
-    graph_attr(crsr.net) <- list()
-    crsr.net <- set_graph_attr(crsr.net, "frag_type", "connection and similarity relations")
+    igraph::E(crsr.net)$type_relation <- NA
+    igraph::E(crsr.net)[ which(crsr.net.edgelist %in% cr.net.edgelist) ]$type_relation <- "cr"
+    igraph::E(crsr.net)[ is.na( igraph::E(crsr.net)$type_relation ) ]$type_relation <- "sr"
+    igraph::graph_attr(crsr.net) <- list()
+    crsr.net <- igraph::set_graph_attr(crsr.net, "frag_type", "connection and similarity relations")
     
     return(crsr.net) 
   }
@@ -198,6 +203,6 @@ make_frag_object <- function(cr, sr, fragments)
     frag_type <- "sr"
   }   
   
-  new(Class="Frag.object", cr.df=cr, sr.df=sr, fragments.df=fragments, frag_type=frag_type)
+  new(Class="Frag.object", df.cr=cr, df.sr=sr, fragments.df=fragments, frag_type=frag_type)
 }
 
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..7c55efb
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,19 @@
+
+
+.check.frag.graph <- function(graph){
+  if(! igraph::is_igraph(graph)) stop("Not a graph object")
+  if(igraph::is_directed(graph)) stop("The 'graph' parameter requires an undirected igraph object.")
+  if(! igraph::is_simple(graph)) stop("Loops and multiple edges are not allowed.")
+}
+
+.check.layer.argument <- function(graph, layer.attr){
+  if(is.null(igraph::vertex_attr(graph, layer.attr)))   stop("'layer.attr' is missing or does not correspond to a vertex attribute of the graph.")
+  if(is.null(igraph::vertex_attr(graph, layer.attr)))   stop("The parameter 'layer.attr' is required.")
+  if (is.null(layer.attr))  stop("No 'layer.attr' argument")
+  if( ! is.character(layer.attr))  stop("The parameter 'layer.attr' requires a character value.")
+  if( ! layer.attr %in% names(igraph::vertex_attr(graph)) ){
+    stop(paste("There is no '", layer.attr, "' vertices attribute.", sep=""))
+  }
+}
+
+
diff --git a/README.Rmd b/README.Rmd
index 678d75e..1b99a07 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -16,7 +16,9 @@ knitr::opts_chunk$set(
 ```
 
 # Archeofrag 
-an R package for refitting and spatial analysis in archeology
+An R package for refitting and spatial analysis in archeology.
+Archeofrag includes methods to analyse fragmented objects in archaeology using refitting relationships between fragments scattered in archaeological spatial units (e.g. stratigraphic layers). Graphs and graph theory are used to model archaeological observations. The package is mainly based on the 'igraph' package for graph analysis. Functions can: 1) create, manipulate, and simulate fragmentation graphs, 2) measure the cohesion and admixture of archaeological spatial units, and 3) characterise the topology of a specific set of refitting relationships. An empirical dataset is also provided as an example.
+
 
 [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
 [![Lifecycle: stable](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://www.tidyverse.org/lifecycle/#maturing)
@@ -43,15 +45,20 @@ an R package for refitting and spatial analysis in archeology
  
 
 # Installation
-The development version can be installed from GitHub with:
+The package can be installed from CRAN with:
+
+```r
+install.packages("archeofrag")
+```
 
+The development version is available on GitHub and can be installed with:
 
 ```r
 # install.packages("devtools")
 devtools::install_github("sebastien-plutniak/archeofrag")
 ```
 
-Note that *Archeofrag* requires the *RBGL* package available through *Bioconductor*:
+Some optional functionalities of *Archeofrag* requires the *RBGL* package available through *Bioconductor*:
 
 ```r
 if (!requireNamespace("BiocManager", quietly = TRUE))
@@ -225,7 +232,7 @@ The `aggreg.factor` parameter affects the distribution of the sizes of the compo
 
 By default, fragments from two spatial units can be disturbed and moved to another other spatial unit. However, the `asymmetric.transport.from` can be used to move fragments from only one given spatial unit.
 
-Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing).
+Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing). Note that the use of this argument requires the RBGL package to be installed.
 
 An example of a complete configuration of the function is:
 
@@ -238,7 +245,7 @@ frag.simul.process(initial.layers=1,
                    components.balance=.4,
                    disturbance=.1,
                    aggreg.factor=0,
-                   planar=T,
+                   planar=TRUE,
                    asymmetric.transport.from="1")
 ```
 
@@ -261,6 +268,7 @@ Setting the  simulator is made easier by using the  `frag.get.parameters` functi
 
 ```{r params}
 params <- frag.get.parameters(abu.g12, layer.attr="layer")
+params
 ```
 
 ```{r simulator-test}
diff --git a/README.md b/README.md
index cd6e123..dc30d47 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,17 @@
+---
+output: github_document
+editor_options: 
+chunk_output_type: console
+---
+
+
+
+
 
 # Archeofrag 
-an R package for refitting and spatial analysis in archeology
+An R package for refitting and spatial analysis in archeology.
+Archeofrag includes methods to analyse fragmented objects in archaeology using refitting relationships between fragments scattered in archaeological spatial units (e.g. stratigraphic layers). Graphs and graph theory are used to model archaeological observations. The package is mainly based on the 'igraph' package for graph analysis. Functions can: 1) create, manipulate, and simulate fragmentation graphs, 2) measure the cohesion and admixture of archaeological spatial units, and 3) characterise the topology of a specific set of refitting relationships. An empirical dataset is also provided as an example.
+
 
 [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
 [![Lifecycle: stable](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://www.tidyverse.org/lifecycle/#maturing)
@@ -27,15 +38,20 @@ an R package for refitting and spatial analysis in archeology
  
 
 # Installation
-The development version can be installed from GitHub with:
+The package can be installed from CRAN with:
+
+```r
+install.packages("archeofrag")
+```
 
+The development version is available on GitHub and can be installed with:
 
 ```r
 # install.packages("devtools")
 devtools::install_github("sebastien-plutniak/archeofrag")
 ```
 
-Note that *Archeofrag* requires the *RBGL* package available through *Bioconductor*:
+Some optional functionalities of *Archeofrag* requires the *RBGL* package available through *Bioconductor*:
 
 ```r
 if (!requireNamespace("BiocManager", quietly = TRUE))
@@ -132,11 +148,11 @@ The `frag.get.layers.pair` function has additional parameters to set the minimum
 ```r
 frag.get.layers.pair(abu.g, layer.attr="layer", sel.layers=c("1", "2"),
                      size.mini=2, mixed.components.only=TRUE)
-#> IGRAPH dd872e8 UN-- 19 22 -- 
+#> IGRAPH 1082003 UN-- 19 22 -- 
 #> + attr: frag_type (g/c), name (v/c), layer (v/c), zmin (v/n), zmax
 #> | (v/n), square (v/c), sherd.type (v/c), thickness (v/n), length (v/n),
 #> | membership (v/n), type_relation (e/c)
-#> + edges from dd872e8 (vertex names):
+#> + edges from 1082003 (vertex names):
 #>  [1] 187--188 165--195 195--196 195--197 196--198 195--204 196--204 197--204
 #>  [9] 198--204 195--25  188--250 27 --28  27 --366 27 --367 28 --367 366--367
 #> [17] 27 --371 332--371 366--371 25 --8   28 --835 835--836
@@ -148,11 +164,11 @@ Additionally, the `frag.get.layers` function can  extract a set of specified spa
 ```r
 frag.get.layers(abu.g, layer.attr="layer", sel.layers="1")
 #> $`1`
-#> IGRAPH 0167b51 UN-- 23 18 -- 
+#> IGRAPH 0cc84c3 UN-- 23 18 -- 
 #> + attr: frag_type (g/c), name (v/c), layer (v/c), zmin (v/n), zmax
 #> | (v/n), square (v/c), sherd.type (v/c), thickness (v/n), length (v/n),
 #> | type_relation (e/c)
-#> + edges from 0167b51 (vertex names):
+#> + edges from 0cc84c3 (vertex names):
 #>  [1] 123--124  187--188  195--196  195--197  196--198  195--204  196--204 
 #>  [8] 197--204  198--204  195--25   301--302  313--314  392--408  435--441 
 #> [15] 477--478  25 --8    435--9999 441--9999
@@ -256,7 +272,7 @@ The `aggreg.factor` parameter affects the distribution of the sizes of the compo
 
 By default, fragments from two spatial units can be disturbed and moved to another other spatial unit. However, the `asymmetric.transport.from` can be used to move fragments from only one given spatial unit.
 
-Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing).
+Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing). Note that the use of this argument requires the RBGL package to be installed.
 
 An example of a complete configuration of the function is:
 
@@ -270,7 +286,7 @@ frag.simul.process(initial.layers=1,
                    components.balance=.4,
                    disturbance=.1,
                    aggreg.factor=0,
-                   planar=T,
+                   planar=TRUE,
                    asymmetric.transport.from="1")
 ```
 
@@ -295,6 +311,30 @@ Setting the  simulator is made easier by using the  `frag.get.parameters` functi
 
 ```r
 params <- frag.get.parameters(abu.g12, layer.attr="layer")
+params
+#> $n.components
+#> [1] 28
+#> 
+#> $vertices
+#> [1] 72
+#> 
+#> $edges
+#> [1] 52
+#> 
+#> $balance
+#> [1] 0.32
+#> 
+#> $components.balance
+#> [1] 0.29
+#> 
+#> $disturbance
+#> [1] 0.04
+#> 
+#> $aggreg.factor
+#> [1] 0.7
+#> 
+#> $planar
+#> [1] FALSE
 ```
 
 
@@ -373,13 +413,13 @@ The `frag.simul.compare` function takes an observed fragmentation graph, generat
 compare.res <- frag.simul.compare(abu.g12, layer.attr="layer",
                                   iter=30, summarise=FALSE)
 head(compare.res$h1.data)
-#>   edges weightsum   balance disturbance   admixture cohesion1 cohesion2
-#> 1    59  344.8793 0.3194444  0.06779661 0.011703301 0.1267001 0.8615966
-#> 2    50  187.0585 0.2916667  0.08000000 0.018833056 0.2100268 0.7711401
-#> 3    55  240.1008 0.3055556  0.14545455 0.057077728 0.1533132 0.7896091
-#> 4    55  268.7331 0.2916667  0.05454545 0.007826262 0.1982092 0.7939646
-#> 5    53  189.6573 0.3055556  0.11320755 0.042060269 0.4362644 0.5216753
-#> 6    56  271.7916 0.3611111  0.08928571 0.021266140 0.2205412 0.7581926
+#>   edges weightsum   balance disturbance  admixture cohesion1 cohesion2
+#> 1    51  173.6126 0.2777778  0.07843137 0.01986406 0.1340509 0.8460851
+#> 2    50  101.8502 0.4444444  0.20000000 0.13116497 0.4130365 0.4557985
+#> 3    50  118.9266 0.1388889  0.14000000 0.05598172 0.1547508 0.7892674
+#> 4    54  176.4861 0.4166667  0.11111111 0.04074675 0.4814336 0.4778196
+#> 5    53  131.0783 0.3472222  0.15094340 0.07717442 0.2598889 0.6629367
+#> 6    51  133.0984 0.3750000  0.11764706 0.04936655 0.2913429 0.6592906
 ```
 
 For each of these parameters, the `frag.simul.summarise` function facilitates the comparison between empirical observed values and simulated values generated for H1 and H2.
@@ -391,13 +431,13 @@ frag.simul.summarise(abu.g12, layer.attr="layer",
                      compare.res$h1.data,
                      compare.res$h2.data)
 #>             H1 != H2? p.value Obs. value/H1 Obs. value/H2
-#> edges           FALSE    0.12         lower        within
-#> weightsum        TRUE    0.02         lower        within
-#> balance         FALSE    0.23        within        within
-#> disturbance     FALSE     0.9         lower         lower
-#> admixture       FALSE    0.97         lower         lower
-#> cohesion1        TRUE       0        higher        within
-#> cohesion2        TRUE       0         lower        within
+#> edges           FALSE    0.67        within        within
+#> weightsum       FALSE    0.11        higher        higher
+#> balance          TRUE       0         lower         lower
+#> disturbance      TRUE       0         lower         lower
+#> admixture        TRUE       0         lower         lower
+#> cohesion1        TRUE       0        higher        higher
+#> cohesion2       FALSE    0.06         lower         lower
 ```
 
 This function returns a data frame with four columns, containing, for each parameter studied:
@@ -503,8 +543,8 @@ rbind(
   "unit1" = frag.cycles(simul.g1, kmax=5),
   "unit2" = frag.cycles(simul.g2, kmax=5))
 #>       3-cycles 4-cycles 5-cycles
-#> unit1       12       10        6
-#> unit2       15        4        1
+#> unit1        0        0        2
+#> unit2        0        1        0
 ```
 
  
@@ -515,11 +555,12 @@ If the `cumulative` parameter is set to `TRUE`, the function returns the cumulat
 
 ```r
 frag.path.lengths(simul.g1)
-#> [1] 31 13  1
+#> [1] 29 50 52 53 30  9
 frag.path.lengths(simul.g2)
-#> [1] 42 18  2
+#> [1]  39  71 102  83  62  32  17   5   1
 frag.path.lengths(simul.g2, cumulative=T)
-#> [1] 1.00000000 0.42857143 0.04761905
+#> [1] 0.382352941 0.696078431 1.000000000 0.813725490 0.607843137 0.313725490
+#> [7] 0.166666667 0.049019608 0.009803922
 ```
 
 In a graph, the shortest path between two vertices is the path including the least number of edges. The diameter of a graph is its longest shortest path.
@@ -528,11 +569,11 @@ The `frag.diameters` function calculates the diameter of each component of the g
 
 ```r
 frag.diameters(simul.g1)
-#> 1 2 3 
-#> 5 4 1
+#> 1 2 3 4 5 6 
+#> 3 0 1 0 0 1
 frag.diameters(simul.g2)
-#> 1 2 3 
-#> 4 5 1
+#> 1 2 3 4 5 6 7 8 9 
+#> 3 1 0 0 0 0 0 0 1
 ```
 
 
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..c765cfe
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,16 @@
+citHeader("To cite archeofrag in publications use:")
+
+citEntry(
+   entry="article",
+   author = "Plutniak, Sébastien",
+   title = "The Strength of Parthood Ties. Modelling Spatial Units and Fragmented Objects with the TSAR Method -- Topological Study of Archaeological Refitting",
+   journal = "Journal of Archaeological Science",
+   year = 2021,
+   doi = "10.1016/j.jas.2021.105501",
+   volume       = "136",
+   number       = "105501",
+   textVersion = "Plutniak, Sébastien (2021). The Strength of Parthood Ties. Modelling Spatial Units and Fragmented Objects with the TSAR Method - Topological Study of Archaeological Refitting. Journal of Archaeological Science. Vol. 136, 105501, DOI: 10.1016/j.jas.2021.105501",
+   header="For the TSAR method:\n"
+)
+
+citFooter("This free open-source software implements academic research by the authors and co-workers. If you use it, please support the project by citing the appropriate journal articles.")
diff --git a/man/Frag.object-class.Rd b/man/Frag.object-class.Rd
index 2eb0ce3..80c1ef6 100644
--- a/man/Frag.object-class.Rd
+++ b/man/Frag.object-class.Rd
@@ -20,8 +20,8 @@
 }
 \section{Slots}{
   \describe{
-    \item{\code{cr.df}:}{Object of class \code{"matrix"} (\code{"data.frame"} are allowed and automatically converted)}
-    \item{\code{sr.df}:}{Object of class \code{"matrix"}  (\code{"data.frame"} are allowed and automatically converted) }
+    \item{\code{df.cr}:}{Object of class \code{"matrix"} (\code{"data.frame"} are allowed and automatically converted)}
+    \item{\code{df.sr}:}{Object of class \code{"matrix"}  (\code{"data.frame"} are allowed and automatically converted)}
     \item{\code{fragments.df}:}{Object of class \code{"data.frame"}}
     \item{\code{frag_type}:}{Object of class \code{"character"} }
   }
@@ -31,14 +31,12 @@
     \item{make_cr_graph}{\code{signature(object = "Frag.object")}:  Makes an undirected graph representing the "connection" relationships between archaeological fragments.}
     \item{make_sr_graph}{\code{signature(object = "Frag.object")}:  Makes an undirected graph representing the "similarity" relationships between archaeological fragments.}
     \item{make_crsr_graph}{\code{signature(object = "Frag.object")}: Makes an undirected graph combining the "connection" and "similarity" relations between archaeological fragments.}
-    \item{show}{\code{signature(object = "Frag.object")}:   \code{show} method for Frag.object }
+    \item{show}{\code{signature(object = "Frag.object")}: \code{show} method for Frag.object }
   }
 }
 
 \author{Sebastien Plutniak <sebastien.plutniak at posteo.net>}
 
-%% ~Make other sections like Warning with \section{Warning }{....} ~
- 
 \seealso{
     \code{\link[archeofrag]{make_frag_object}},
     \code{\link[archeofrag]{make_cr_graph}},
diff --git a/man/LiangAbu.Rd b/man/LiangAbu.Rd
index 721608d..74347f4 100644
--- a/man/LiangAbu.Rd
+++ b/man/LiangAbu.Rd
@@ -14,18 +14,18 @@
         \item A similarity relationship between fragments is defined if there is an acceptable likelihood that those fragments were part of the same object.
     }
     
-    The dataset is composed of three tables, \code{cr.df}, \code{sr.df}, \code{fragments.info}.
+    The dataset is composed of three tables, \code{df.cr}, \code{df.sr}, \code{fragments.info}.
     \itemize{
-        \item\code{cr.df}: "connection" relationships between fragments.
-        \item\code{sr.df}: "similarity" relationships between fragments.
+        \item\code{df.cr}: "connection" relationships between fragments.
+        \item\code{df.sr}: "similarity" relationships between fragments.
         \item\code{fragments.info}: contextual information concerning each fragment.
     }
 }
 \usage{data(LiangAbu)}
 \format{
     \itemize{
-        \item\code{cr.df} is a 56x2 matrix. Each line describes a connection relationship between two fragments. There respective unique identifiers are given in column "frg_id1" and in column "frg_id2". 
-        \item\code{sr.df} is a 147x2 matrix. Column "frg_id" gives a fragment unique identifier, column "su_id" gives a unique identifier for the group of similar fragments it belongs to (similarity unit).
+        \item\code{df.cr} is a 56x2 matrix. Each line describes a connection relationship between two fragments. There respective unique identifiers are given in column "frg_id1" and in column "frg_id2". 
+        \item\code{df.sr} is a 147x2 matrix. Column "frg_id" gives a fragment unique identifier, column "su_id" gives a unique identifier for the group of similar fragments it belongs to (similarity unit).
         \item\code{fragments.info} is 177x8 data frame:
         \itemize{
             \item frg_id: unique fragment identifier 
diff --git a/man/figures/README-manipulate-plot-abu-1.png b/man/figures/README-manipulate-plot-abu-1.png
index fa11dcf..74de551 100644
Binary files a/man/figures/README-manipulate-plot-abu-1.png and b/man/figures/README-manipulate-plot-abu-1.png differ
diff --git a/man/figures/README-manipulate-plot-abu2-1.png b/man/figures/README-manipulate-plot-abu2-1.png
index f4894dc..0af832b 100644
Binary files a/man/figures/README-manipulate-plot-abu2-1.png and b/man/figures/README-manipulate-plot-abu2-1.png differ
diff --git a/man/figures/README-simulator-test2-admix-1.png b/man/figures/README-simulator-test2-admix-1.png
index 563ee38..995dd01 100644
Binary files a/man/figures/README-simulator-test2-admix-1.png and b/man/figures/README-simulator-test2-admix-1.png differ
diff --git a/man/figures/README-simulator-test2-edges-1.png b/man/figures/README-simulator-test2-edges-1.png
index aa72b24..f5dd4e8 100644
Binary files a/man/figures/README-simulator-test2-edges-1.png and b/man/figures/README-simulator-test2-edges-1.png differ
diff --git a/man/frag.edges.weighting.Rd b/man/frag.edges.weighting.Rd
index 2c153ff..94476b2 100644
--- a/man/frag.edges.weighting.Rd
+++ b/man/frag.edges.weighting.Rd
@@ -56,6 +56,7 @@ The graph, with an additional "weight" edge attribute and, if the distance has b
 g <- frag.simul.process(n.components=20, vertices=50, disturbance=.15)
 frag.edges.weighting(g , "layer")
 # with morphometric and spatial parameters:
+library(igraph)
 V(g)$morpho <- sample(1:20, 50, replace=TRUE)
 V(g)$x <- sample(1:100, 50, replace=TRUE)
 V(g)$y <- sample(1:100, 50, replace=TRUE)
diff --git a/man/frag.get.layers.Rd b/man/frag.get.layers.Rd
index cee4cc1..473f8ed 100644
--- a/man/frag.get.layers.Rd
+++ b/man/frag.get.layers.Rd
@@ -25,7 +25,7 @@ A list with a graph for each selected stratigraphic layer.
 
 \examples{
 g <- frag.simul.process(n.components=20, vertices=50, disturbance = .15)
-V(g)$layers <- c(rep("layer1", 20), rep("layer2", 20), rep("layer3", 10))
+igraph::V(g)$layers <- c(rep("layer1", 20), rep("layer2", 20), rep("layer3", 10))
 frag.get.layers(g, layer.attr="layers", sel.layers=c("layer1", "layer2"))
 }
 \keyword{extraction}
diff --git a/man/frag.get.layers.pair.Rd b/man/frag.get.layers.pair.Rd
index d504ac7..618ab44 100644
--- a/man/frag.get.layers.pair.Rd
+++ b/man/frag.get.layers.pair.Rd
@@ -28,7 +28,7 @@ An undirected graph object.
 
 \examples{
 g <- frag.simul.process(n.components=20, vertices=50, disturbance=.15)
-V(g)$layers <- c(rep("layer1", 20), rep("layer2", 20), rep("layer3", 10))
+igraph::V(g)$layers <- c(rep("layer1", 20), rep("layer2", 20), rep("layer3", 10))
 
 frag.get.layers.pair(g, "layers", sel.layers=c("layer2","layer3"),
                      size.mini=2, mixed.components.only=FALSE)
diff --git a/man/frag.get.parameters.Rd b/man/frag.get.parameters.Rd
index c6ef624..b59cecf 100644
--- a/man/frag.get.parameters.Rd
+++ b/man/frag.get.parameters.Rd
@@ -13,13 +13,19 @@
 This function is a convenient function to obtain general information about a fragmentation graph. It is particularly useful for setting the parameters of the \code{frag.simul.process} function.
 It returns the number of components, vertices, and edges, the balance (proportion of fragments in the smaller layer), components balance (proportion of components in the poorest layer), the disturbance, the aggregation factor, and if the graph is planar or not.
 The aggregation factor reflects the diversity of the components' edge counts. The factor is calculated by: 1 - 1/(1 + sd(edge counts of the components)).
+The optional RBGL package is required to determine the planarity of the graph. If it is not installed, the `planar` value is set to FALSE by default.
 }
 \value{
 A list with the values for the parameters.
 }
 \author{Sebastien Plutniak <sebastien.plutniak at posteo.net>}
 
-\seealso{\link{frag.get.layers.pair}, \link{frag.simul.process}}
+\seealso{
+		 \link[archeofrag]{frag.get.layers.pair}, 
+		 \link[archeofrag]{frag.simul.process}, 
+		 \link[stats]{sd}, 
+		 \code{\link[RBGL]{boyerMyrvoldPlanarityTest}}
+		 }
 
 \examples{
 g <- frag.simul.process(n.components=20, vertices=50)
diff --git a/man/frag.simul.compare.Rd b/man/frag.simul.compare.Rd
index 715685c..c6ab25a 100644
--- a/man/frag.simul.compare.Rd
+++ b/man/frag.simul.compare.Rd
@@ -3,13 +3,14 @@
 
 \title{From an observed fragmentation graph, simulates two series of graphs corresponding to two deposition hypotheses.}
 \description{Given an observed fragmentation graph, simulates two series of graphs corresponding to two deposition hypotheses, compares their properties and returns a summary table.}
-\usage{frag.simul.compare(graph, layer.attr, iter, summarise=TRUE)}
+\usage{frag.simul.compare(graph, layer.attr, iter, summarise=TRUE, ...)}
 
 \arguments{
   \item{graph}{An undirected \code{igraph} object. The 'observed' graph to compare to simulated graphs.} 
   \item{layer.attr}{Character. The name of the vertices attribute giving the layer of the fragments.}
   \item{iter}{Numerical. The number of simulated graphs to generate for each hypothesis (minimal value: 30).}
   \item{summarise}{Logical. Whether to report a comparative summary of the results.}
+  \item{...}{Further arguments passed to the `frag.simul.process` function.}
 }
 \details{
 This function is a convenient wrapper integrating several functions of the \code{archeofrag} package to compare an observed fragmentation graph to similar simulated graphs. The \code{frag.simul.process} is used to generate two series of graphs from the properties of the observed graph: the first series is generated under the formation hypothesis H1 (one initial spatial unit) and the second series is generated under the hypothesis H2 (two initial spatial units). 
@@ -30,7 +31,7 @@ A named list with three items:
 
 \seealso{
 \code{\link[archeofrag]{frag.simul.process}},
-\code{\link[archeofrag]{frag.simul.summarise}},
+\code{\link[archeofrag]{frag.simul.summarise}}
 }
 \examples{
 g <- frag.simul.process(n.components=20, vertices=50, disturbance=.15)
diff --git a/man/frag.simul.process.Rd b/man/frag.simul.process.Rd
index cea4c6a..2baab43 100644
--- a/man/frag.simul.process.Rd
+++ b/man/frag.simul.process.Rd
@@ -9,7 +9,7 @@ Simulate the fragmentation of archaeological objects scattered in two stratigrap
 \usage{
 frag.simul.process(initial.layers=2, n.components, vertices=Inf,
                    edges=Inf, balance=.5, components.balance=.5,
-                   disturbance=0, aggreg.factor=0, planar=TRUE,
+                   disturbance=0, aggreg.factor=0, planar=FALSE,
                    asymmetric.transport.from=NULL, 
                    from.observed.graph=NULL, observed.layer.attr=NULL)
 }
@@ -34,6 +34,8 @@ Fragments are represented by vertices and the "connection" relationships ("refit
 
 Some parameters are optional or depend on other parameters (messages are displayed accordingly). Namely, if two \code{initial.layers} are set, then only one of the \code{vertices} and \code{edges} parameters can be used. Using only one layer as the initial condition enables to constraint the graph with the number of vertices only, the number of edges only, or both. 
 
+Note that using the \code{planar} argument requires to install the optional RBGL package.
+
 The \code{disturbance} determines the proportion of fragments to "move" from one layer to another. Consequently, it generates inter-layer relationships.
 Note that the \code{balance} parameter determines the proportion of fragments in the first layer before the application of the disturbance process.
 
@@ -48,9 +50,9 @@ An igraph object with a "frag_type" graph attribute (with the value "cr", for "c
 }
 \author{Sebastien Plutniak <sebastien.plutniak at posteo.net>}
 \seealso{
-    \code{\link[RBGL]{boyerMyrvoldPlanarityTest}},
     \code{\link[archeofrag]{frag.get.parameters}},
-    \code{\link[archeofrag]{frag.edges.weighting}}
+    \code{\link[archeofrag]{frag.edges.weighting}},
+    \code{\link[RBGL]{boyerMyrvoldPlanarityTest}}
 }
 
 \examples{
@@ -63,8 +65,8 @@ g <- frag.simul.process(initial.layers=1,
                             balance=.5,
                             components.balance=.5,
                             disturbance=.1,
-                            planar=TRUE)
-plot(g, vertex.color=factor(V(g)$layer), 
+                            planar=FALSE)
+plot(g, vertex.color=factor(igraph::V(g)$layer), 
      vertex.size=4, vertex.label=NA)
 
 }
diff --git a/vignettes/archeofrag-vignette.Rmd b/vignettes/archeofrag-vignette.Rmd
index 2761652..cb05329 100644
--- a/vignettes/archeofrag-vignette.Rmd
+++ b/vignettes/archeofrag-vignette.Rmd
@@ -197,7 +197,7 @@ The `aggreg.factor` parameter affects the distribution of the sizes of the compo
 
 By default, fragments from two spatial units can be disturbed and moved to another other spatial unit. However, the `asymmetric.transport.from` can be used to move fragments from only one given spatial unit.
 
-Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing).
+Finally, the `planar` argument determines if the generated graph has to be planar or not (a graph is planar when it can be drawn on a plane, without edges crossing). This function requires to install the optional `RBGL` package.
 
 An example of a complete configuration of the function is:
 
@@ -210,7 +210,7 @@ frag.simul.process(initial.layers=1,
                    components.balance=.4,
                    disturbance=.1,
                    aggreg.factor=0,
-                   planar=T,
+                   planar=FALSE,
                    asymmetric.transport.from="1")
 ```
 
@@ -233,6 +233,7 @@ Setting the  simulator is made easier by using the  `frag.get.parameters` functi
 
 ```{r params}
 params <- frag.get.parameters(abu.g12, layer.attr="layer")
+params
 ```
 
 ```{r simulator-test}
@@ -282,7 +283,7 @@ par(mar=c(5, 4, 4, 2))
 ```{r simulator-test2-edges, message=FALSE, fig.align="center", fig.width=4, fig.height=3}
 edges.res <- sapply(test2.results,
                     function(g) frag.get.parameters(g, "layer")$edges)
-plot(density(edges.res), main="Edges")
+plot(stats::density(edges.res), main="Edges")
 abline(v=params$edges, col="red")
 ```
 
@@ -291,12 +292,12 @@ Similarly, the empirical admixture value is lower than the simulated admixture v
 ```{r simulator-test2-admix, message=FALSE, fig.align="center", fig.width=4, fig.height=3}
 admix.res <- sapply(test2.results,
                     function(g) frag.layers.admixture(g, "layer"))
-plot(density(admix.res), main="Admixture")
+plot(stats::density(admix.res), main="Admixture")
 abline(v=frag.layers.admixture(abu.g12, "layer"), col="red")
 ```
  
 Two functions (`frag.simul.compare` and `frag.simul.summarise`) facilitate the  execution of the analytical process described above on the initial number of spatial units.
-The `frag.simul.compare` function takes an observed fragmentation graph, generates two series of simulated graphs corresponding to two hypotheses on the number of initial spatial units (H1 for 1 initial spatial unit and H2 for two initial spatial units), and returns a data frame of measurements made on each series (including the edge count, weights sum, balance value, disturbance value, admixture value, and cohesion values of the two spatial units).
+The `frag.simul.compare` function takes an observed fragmentation graph, generates two series of simulated graphs corresponding to two hypotheses on the number of initial spatial units (H1 for one initial spatial unit and H2 for two initial spatial units), and returns a data frame of measurements made on each series (including the edge count, weights sum, balance value, disturbance value, admixture value, and cohesion values of the two spatial units).
 
 
 ```{r simul-compare, message=FALSE}
@@ -305,6 +306,7 @@ compare.res <- frag.simul.compare(abu.g12, layer.attr="layer",
 head(compare.res$h1.data)
 ```
 
+
 For each of these parameters, the `frag.simul.summarise` function facilitates the comparison between empirical observed values and simulated values generated for H1 and H2.
 
 
@@ -315,7 +317,7 @@ frag.simul.summarise(abu.g12, layer.attr="layer",
 ```
 
 This function returns a data frame with four columns, containing, for each parameter studied:
-
+  
 1. whether the series of H1 values are statistically different to the H2 series (Boolean), 
 2. the p-value of the Wilcoxon test (numerical), 
 3. whether the observed value is "within", "higher", or "lower" to the interquartile range of values for H1, 
@@ -339,7 +341,7 @@ abu.sr <- make_sr_graph(abu.frag)
 ```
 
 The `frag.relations.by.layers` function returns a table with the number of similarity relationships in and between spatial units, e.g., in the top three layers at Liang Abu:
-
+  
 ```{r count-similarity}
 # count of similarity relationships in and between layers:
 simil.by.layers.df <- frag.relations.by.layers(abu.sr, "layer")
@@ -347,7 +349,7 @@ simil.by.layers.df
 ```
 
 These values can be observed as percentages:
-
+  
 ```{r similarity-perc-tab}
 # percentage of similarity relationships in and between layers:
 round(simil.by.layers.df / sum(simil.by.layers.df, na.rm=T) * 100, 0)
@@ -361,13 +363,13 @@ The expected result is observed for Liang Abu surface and the first two layers,
 simil.dist <- max(c(simil.by.layers.df), na.rm=T) - simil.by.layers.df
 simil.dist <- as.dist(simil.dist)
 # hierarchical clustering:
-clust.res <- hclust(simil.dist, method="ward.D2")
+clust.res <- stats::hclust(simil.dist, method="ward.D2")
 ```
 
 ```{r similarity-dendr-fig, fig.width = 3, fig.height = 3, fig.align="center", fig.cap="Hierarchical clustering of the pottery layers in Liang Abu (distance: based on the number of similarity relationships; clustering method: Ward).", eval=T, message=F}
 clust.res$labels <- as.character(factor(clust.res$labels, 
-                     levels=c("0", "1", "2"),
-                     labels=c("layer 0", "layer 1", "layer 2")))
+                                        levels=c("0", "1", "2"),
+                                        labels=c("layer 0", "layer 1", "layer 2")))
 plot(clust.res, hang=-1, axes=F, ann=F)
 ```
 
@@ -392,14 +394,14 @@ In a graph, a cycle is a path in which only the first and last vertices are repe
 The `frag.cycles` function searches for cycles in a graph and returns the number of cycles found for different cycle lengths. 
 The `kmax` parameter determines the maximal length of the cycles to search for.
 Let us compare the cycles found in the two spatial units of the artificial graph:
-
+  
 ```{r cycles-simul1}
 rbind(
   "unit1" = frag.cycles(simul.g1, kmax=5),
   "unit2" = frag.cycles(simul.g2, kmax=5))
 ```
 
- 
+
 The `frag.path.lengths` function returns the distribution of the path lengths in the graph (i.e., the number of edges between each pair of vertices). 
 This function returns a vector whose first element is the frequency of the paths of length 1, the second element is the frequency of the paths of length 2, etc. 
 If the `cumulative` parameter is set to `TRUE`, the function returns the cumulative relative frequency of the path lengths.
@@ -419,4 +421,3 @@ frag.diameters(simul.g2)
 ```
 
 
-