diff --git a/R/RN_dims.R b/R/RN_dims.R
new file mode 100644
index 0000000..4ce856a
--- /dev/null
+++ b/R/RN_dims.R
@@ -0,0 +1,73 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' # RN_dims_UNI() modified to enter the interaction and cover data
+#The size of a network can be described by the number of nodes (*N*) and links (*L*), and its complexity is #proportional to network connectance. *C* can be estimated on different ways depending on the type of #network, but it always measures the proportion of links relative to the maximum number of links that could #be possible in the network. In the case of recruitment networks we use the formula $C = L /(N^2-N)$ since #the node "Open" does not act as a recruit (i.e. Open is represented by a row of zeroes in the adjacency #matrix).
+#For facilitation and competition that are bipartite networks, the nodes in eahc guild are #provided( nurse/facilitated, Canopy (depressor)/recruit) seprately, and open is not a node
+#'
+RN_dims <- function(int_data,cover_data, int_type=c("rec","fac","comp")){
+
+  if(int_type=="rec"){
+
+    int_data <-comm_to_RN_UNI(int_data,cover_data)
+
+    # FUNCTION
+    df <- int_data
+    n_nodes <- length(unique(c(df$Canopy, df$Recruit)))
+    n_links <- sum(df$Pcr)
+    connectance <- n_links/(n_nodes^2 - n_nodes)
+
+    out <- data.frame(c(n_nodes, n_links, connectance))
+    colnames(out) <- c("Value")
+    rownames(out) <- c("Num. Nodes", "Num. Links", "Connectance")
+    return(out)
+
+  }
+
+  if(int_type=="fac"){
+
+    int_data <-suppressWarnings(RN_to_matrix(int_data,cover_data, int_type="fac",weight="Pcr"))
+
+    df <-int_data
+    n_nodes_nurse <-dim(df)[2]
+    n_nodes_facilitated <-dim(df)[1]
+    n_links <- sum(df)
+    connectance <-n_links/(n_nodes_nurse*n_nodes_facilitated)
+
+    out <- data.frame(c(n_nodes_nurse,n_nodes_facilitated, n_links, connectance))
+    colnames(out) <- c("Value")
+    rownames(out) <- c("Num. Nurse sp","Num. Facilitated sp", "Num. Links", "Connectance")
+    return(out)
+
+  }
+
+
+  if(int_type=="comp"){
+
+    int_data <-suppressWarnings(RN_to_matrix(int_data,cover_data, int_type="comp",weight="Pcr"))
+
+    df <-int_data
+    n_nodes_canopy_depressing <-dim(df)[2]
+    n_nodes_recruit_depressed <-dim(df)[1]
+    n_links <- sum(df)
+    connectance <-n_links/(n_nodes_canopy_depressing*n_nodes_recruit_depressed)
+
+    out <- data.frame(c(n_nodes_canopy_depressing, n_nodes_recruit_depressed, n_links, connectance))
+    colnames(out) <- c("Value")
+    rownames(out) <- c("Num. Canopy depressing sp","Num. Recruit depressed sp", "Num. Links", "Connectance")
+    return(out)
+
+
+  }
+
+
+}
+
diff --git a/R/RN_heatmap.R b/R/RN_heatmap.R
new file mode 100644
index 0000000..5d4515a
--- /dev/null
+++ b/R/RN_heatmap.R
@@ -0,0 +1,174 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#' @param weight
+#' @param scale_top
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+RN_heatmap <- function(int_data,cover_data,int_type=c("rec","fac","comp"), weight = c("Pcr","Fcr","Dcr","Dro","Ns", "NintC", "NintA", "RII"), scale_top = 1) {
+
+  if(int_type=="rec"){
+
+    index<-suppressWarnings(associndex(int_data,cover_data,expand="yes",rm_sp_no_cover="allsp"))
+    data<-comm_to_RN_UNI(int_data,cover_data)
+    data$Dcr<-data$Fcr/data$Ac
+
+    # calculate Dro for each recuit species
+    open_rows <- data[data$Canopy == "Open", ]
+    ratios <- tapply(open_rows$Fcr / open_rows$Ac, open_rows$Recruit, mean)
+    data$Dro <- ratios[data$Recruit]
+
+
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    index$int<-paste(index$Canopy,index$Recruit, sep="_")
+    db<-merge(data[,c("int","Canopy","Recruit","Fcr","Icr","Pcr","Dcr","Dro")],index[,c("int","Ns", "NintC", "NintA", "RII")],   by="int", all.x=T)
+    db[is.na(db)]<-0
+
+    int_data<-db
+
+    if (weight %in% c("Ns", "NintC", "NintA", "RII")) {
+      warning("Since the index specified in the 'weight' argument is defined relative to 'Open', rows or columns labeled 'Open' are mathematically zero and not biologically meaningful")}
+
+    require(ggplot2)
+    # manually set node order
+    canopy_order <- unique(int_data$Canopy)
+    canopy_order <- canopy_order[!canopy_order %in% c('Open')]
+    canopy_order <- c("Open", canopy_order)
+    int_data$Canopy2 <- factor(int_data$Canopy, levels = canopy_order)
+    recruit_order <- sort(unique(int_data$Canopy), decreasing = TRUE)
+    recruit_order <- recruit_order[!recruit_order %in% c('Open')]
+    recruit_order <- c(recruit_order, "Open")
+    int_data$Recruit2 <- factor(int_data$Recruit, levels = recruit_order)
+
+
+    # Make weight variable
+    int_data$weight_v <- as.numeric(int_data[[weight]])
+
+    # Lowest (non-zero) and highest values of the weighting variable
+    highest_W <- max(int_data$weight_v)
+    lowest_W <- min(int_data$weight_v[int_data$weight_v>0])
+
+    # Plot the heatmap
+    p<-ggplot(int_data, aes(Canopy2, Recruit2, fill= weight_v)) +
+      geom_tile(aes(height = 1, alpha = weight_v != 0),
+                colour = "gray", linewidth = 0.25) +
+      scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0), guide = "none") +
+      scale_x_discrete(position = "top") +
+      labs(fill = weight,x = "Canopy", y = "Recruit") +
+      theme(panel.grid = element_blank(),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0))
+
+    return(p)
+
+  }
+
+
+
+  if(int_type=="fac"){
+
+    data<-comm_to_RN_BI(int_data,cover_data)
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    df<-suppressWarnings(int_significance_BI (data))
+    df<-df[df$Effect_int=="Enhancing",]
+    fac_int<-paste(df$Canopy,df$Recruit, sep="_")
+    db<-suppressWarnings(associndex_UNISITE_BI(comm_to_RN_BI(int_data,cover_data)))
+    db$int<-paste(db$Canopy,db$Recruit, sep="_")
+    db<-merge(db, data[,c("int","Icr","Pcr")], by="int", all.x=T)
+    db<-db[db$int%in%fac_int,]
+
+    int_data<-db
+
+    require(ggplot2)
+    # manually set node order
+    canopy_order <- unique(int_data$Canopy)
+    int_data$Canopy2 <- factor(int_data$Canopy, levels = canopy_order)
+    recruit_order <- sort(unique(int_data$Recruit), decreasing = TRUE)
+    int_data$Recruit2 <- factor(int_data$Recruit, levels = recruit_order)
+
+    # Add recruitment density as another weighting variable
+    int_data$Dcr <- int_data$Fcr/int_data$Ac
+
+    # Make weight variable
+    int_data$weight_v <- as.numeric(int_data[[weight]])
+
+    # Lowest (non-zero) and highest values of the weighting variable
+    highest_W <- max(int_data$weight_v)
+    lowest_W <- min(int_data$weight_v[int_data$weight_v>0])
+
+    # Plot the heatmap
+    p<-ggplot(int_data, aes(Canopy2, Recruit2, fill= weight_v)) +
+      geom_tile(aes(height = 1, alpha = weight_v != 0),
+                colour = "gray", linewidth = 0.25) +
+      scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0), guide = "none") +
+      scale_x_discrete(position = "top") +
+      labs(fill = weight,x = "Canopy", y = "Recruit") +
+      theme(panel.grid = element_blank(),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0))
+
+    return(p)
+
+  }
+
+
+  if(int_type=="comp"){
+
+
+    data<-comm_to_RN_UNI_COMP(int_data,cover_data)
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    df<-suppressWarnings(int_significance_BI_COMP (data))
+    df<-df[df$Effect_int=="Depressing",]
+    comp_int<-paste(df$Canopy,df$Recruit, sep="_")
+    db<-suppressWarnings(associndex_UNISITE_BI_COMP(comm_to_RN_UNI_COMP(int_data,cover_data)))
+    db$int<-paste(db$Canopy,db$Recruit, sep="_")
+    db<-merge(db, data[,c("int","Icr","Pcr")], by="int", all.x=T)
+    db<-db[db$int%in%comp_int,]
+
+    int_data<-db
+
+    require(ggplot2)
+    # manually set node order
+    canopy_order <- unique(int_data$Canopy)
+    int_data$Canopy2 <- factor(int_data$Canopy, levels = canopy_order)
+    recruit_order <- sort(unique(int_data$Recruit), decreasing = TRUE)
+    int_data$Recruit2 <- factor(int_data$Recruit, levels = recruit_order)
+
+    # Add recruitment density as another weighting variable
+    int_data$Dcr <- int_data$Fcr/int_data$Ac
+
+    # Make weight variable
+    int_data$weight_v <- as.numeric(int_data[[weight]])
+
+    if(weight%in% c("Fcr","Icr","Pcr","Dcr","Dro")){
+      # Lowest (non-zero) and highest values of the weighting variable
+      highest_W <- max(int_data$weight_v)
+      lowest_W <- min(int_data$weight_v[int_data$weight_v>0])
+    }
+
+
+    if(weight %in% c("Ns", "NintC", "NintA", "RII")){
+      nonzero_vals <- int_data$weight_v[int_data$weight_v != 0]
+      highest_W <- max(nonzero_vals)
+      lowest_W <- min(nonzero_vals)
+    }
+    # Plot the heatmap
+    p<-ggplot(int_data, aes(Canopy2, Recruit2, fill= weight_v)) +
+      geom_tile(aes(height = 1, alpha = weight_v != 0),
+                colour = "gray", linewidth = 0.25) +
+      scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0), guide = "none") +
+      scale_x_discrete(position = "top") +
+      labs(fill = weight,x = "Canopy", y = "Recruit") +
+      theme(panel.grid = element_blank(),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0))
+
+    return(p)
+
+
+  }
+
+}
+
+
diff --git a/R/RN_to_matrix.R b/R/RN_to_matrix.R
new file mode 100644
index 0000000..018fa54
--- /dev/null
+++ b/R/RN_to_matrix.R
@@ -0,0 +1,87 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#' @param weight
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#' #This function returns a matrix with canopy species as columns and recruits as rows, containing the association index specified in the weight argument. If int_type is set to "rec," all possible interactions between species with known coverage are displayed. If set to "fac" (or "comp"), a nonzero value appears when recruits associate with canopy plants more (or less) than expected based on their coverage; otherwise, the value is zero.
+#'
+RN_to_matrix<- function(int_data,cover_data, int_type=c("rec", "fac","comp"), weight = c("Fcr","Dcr","Dro","Ns", "NintC", "NintA", "RII")){
+
+  if (!"Open" %in% int_data$Canopy) stop("ERROR: your data does not contain a node named Open or it is spelled differently.")
+
+  if(int_type=="rec"){
+
+    index<-suppressWarnings(associndex(int_data,cover_data,expand="yes",rm_sp_no_cover="allsp"))
+    data<-comm_to_RN_UNI(int_data,cover_data)
+    data$Dcr<-data$Fcr/data$Ac
+
+    # calculate Dro for each recruit species
+    open_rows <- data[data$Canopy == "Open", ]
+    ratios <- tapply(open_rows$Fcr / open_rows$Ac, open_rows$Recruit, mean)
+    data$Dro <- ratios[data$Recruit]
+
+
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    index$int<-paste(index$Canopy,index$Recruit, sep="_")
+    db<-merge(data[,c("int","Canopy","Recruit","Fcr","Icr","Pcr","Dcr","Dro")],index[,c("int","Ns", "NintC", "NintA", "RII")],   by="int", all.x=T)
+    db[is.na(db)]<-0
+    net<-suppressWarnings(RN_to_matrix_UNI(db, weight))
+
+    if (weight %in% c("Ns", "NintC", "NintA", "RII")) {
+      warning("Since the index specified in the 'weight' argument is defined relative to 'Open', rows or columns labeled 'Open' are mathematically zero and not biologically meaningful")
+    }
+  }
+
+  if(int_type=="fac"){
+
+    data<-comm_to_RN_BI(int_data,cover_data)
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    df<-suppressWarnings(int_significance_BI (data))
+
+    if(nrow(df[df$Effect_int=="Enhancing",])==0){
+
+      stop("There is not any recruitment enhancing interaction")
+    }else{
+
+      df<-df[df$Effect_int=="Enhancing",]
+      fac_int<-paste(df$Canopy,df$Recruit, sep="_")
+      db<-suppressWarnings(associndex_UNISITE_BI(comm_to_RN_BI(int_data,cover_data)))
+      db$int<-paste(db$Canopy,db$Recruit, sep="_")
+      db<-merge(db, data[,c("int","Icr","Pcr")], by="int", all.x=T)
+      db<-db[db$int%in%fac_int,]
+      net<-suppressWarnings(RN_to_matrix_BI(db, weight))
+    }
+
+  }
+
+  if(int_type=="comp"){
+
+    data<-comm_to_RN_UNI_COMP(int_data,cover_data)
+    data$int<-paste(data$Canopy,data$Recruit, sep="_")
+    df<-suppressWarnings(int_significance_BI_COMP (data))
+
+    if(nrow(df[df$Effect_int=="Depressing",])==0){
+
+      stop("There is not any recruitment depressing interaction")
+    }else{
+
+      df<-df[df$Effect_int=="Depressing",]
+      comp_int<-paste(df$Canopy,df$Recruit, sep="_")
+      db<-suppressWarnings(associndex_UNISITE_BI_COMP(comm_to_RN_UNI_COMP(int_data,cover_data)))
+      db$int<-paste(db$Canopy,db$Recruit, sep="_")
+      db<-merge(db, data[,c("int","Icr","Pcr")], by="int", all.x=T)
+      db<-db[db$int%in%comp_int,]
+      net<-suppressWarnings(RN_to_matrix_BI(db, weight))
+    }
+  }
+
+
+  return(net)
+}
+
diff --git a/R/associndex.R b/R/associndex.R
new file mode 100644
index 0000000..3c74fca
--- /dev/null
+++ b/R/associndex.R
@@ -0,0 +1,78 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param expand
+#' @param rm_sp_no_cover
+#' @param threshold_density
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+#Calculate different association indices considering different interaction types refer to:
+#expand:expands non observed interactions (yes or no)
+#rm_sp_no_cover: "allsp" removes any species without cover data, while "onlycanopy":removes only canopy species without cover
+#'
+associndex<- function(int_data, cover_data, expand=c("yes","no"), rm_sp_no_cover=c("allsp","onlycanopy"),
+                      threshold_density=NULL)
+
+{
+
+  if (!"Open" %in% int_data$Canopy) stop("ERROR: your data does not contain a node named Open or it is spelled differently. Data for recruitment in Open is required to calculate the indices.")
+
+  # Suggest a threshold value based on recruitment density outliers according to a Weibull distribution. Uses "extremevalues" package.
+  int_data_0<-comm_to_RN_UNI(int_data, cover_data)
+  db_inter_0 <- pre_associndex_UNISITE_UNI(int_data_0)
+  db_inter_0$Dcr <- db_inter_0$Fcr/db_inter_0$Ac
+  db_inter_0$Dro <- db_inter_0$Fro/db_inter_0$Ao
+  y <- rbind(db_inter_0$Dcr,db_inter_0$Dro)
+  y <- y[which(y>0)] # Zero densities are common, so we will check only for suspiciously high density values.
+  K <- getOutliers(y, method = "I", distribution = "weibull", rho=c(1,1)) # We consider that a value is an outlier
+  # if it is above the limit where less
+  # than 1 observation is expected.
+
+  if(is.null(threshold_density)){
+    threshold_density = K$limit[[2]]
+    message("Based on Weibull distribution fitted to the observed values, the threshold density has been set to: ", threshold_density, ".")
+  }
+
+
+  if(expand=="yes" & rm_sp_no_cover=="allsp"){int_type ="rec"}
+  if(expand=="yes" & rm_sp_no_cover=="onlycanopy"){int_type ="comp"}
+  if(expand=="no" & rm_sp_no_cover=="onlycanopy"){int_type ="fac"}
+  if(expand=="no" & rm_sp_no_cover=="allsp"){int_type ="unused"}
+
+  if(int_type=="rec"){
+
+    int_data<-comm_to_RN_UNI(int_data, cover_data)
+    db_inter<-associndex_UNISITE_UNI(int_data, threshold_density)
+
+  }
+
+  if(int_type=="fac"){
+
+    int_data<-comm_to_RN_BI(int_data, cover_data)
+    db_inter<-associndex_UNISITE_BI(int_data, threshold_density)
+
+  }
+
+  if(int_type=="comp"){
+
+    int_data<-comm_to_RN_UNI_COMP(int_data, cover_data)
+    db_inter<-associndex_UNISITE_BI_COMP(int_data, threshold_density)
+
+  }
+
+  if(int_type=="unused"){
+
+    db_inter<-NULL
+    warning("The object is not provided because the combination of arguments is not commonly used in plant-plant interaction networks")
+  }
+
+
+  return(db_inter)
+}
+
diff --git a/R/canopy_service_test.R b/R/canopy_service_test.R
new file mode 100644
index 0000000..65db051
--- /dev/null
+++ b/R/canopy_service_test.R
@@ -0,0 +1,89 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' #como canopy_service_test_UNI pero sin quitar especies de reclutas
+#que no tienen cover ( solo las canopy que no tienen cover)
+#sustituimos pre_associndex_UNISITE_UNI por
+#associndex(int_data,cover_data, expand="yes", rm_sp_no_cover="onlycanopy")
+#calcula para cada especie de nodriza, si alberga bajo su copa una mayor (o menor)
+# densidad de reclutas (de cualquier especie) que la esperada en funcion del area que ocupa.
+#The input data are two files: the field data containing interactions and species cover
+#'
+canopy_service_test <- function(int_data,cover_data){
+
+  df<-associndex(int_data,cover_data, expand="yes", rm_sp_no_cover="onlycanopy")
+  sp_Fc <- stats::aggregate(Fcr ~ Canopy, data = df, FUN = sum)
+  sp_Ac <- stats::aggregate(Ac ~ Canopy, data = df, FUN = max)
+  sp_Fro <- stats::aggregate(Fro ~ Canopy, data = df, FUN = sum)
+  sp_Ao <- stats::aggregate(Ao ~ Canopy, data = df, FUN = max)
+  n_tests <- dim(sp_Fc)[1]
+  df <- data.frame(c(sp_Fc, sp_Ac, sp_Fro, sp_Ao))
+  myvars <- names(df) %in% c("Canopy.1", "Canopy.2", "Canopy.3")
+  df <- df[!myvars]
+  colnames(df) <- c("Canopy", "Fc", "Ac", "Fro", "Ao")
+  df$exp_p <- df$Ac/(df$Ac+df$Ao) # Expected probability of success (i.e. of recruiting under canopy)
+
+  # Testability through Binomial test
+
+  df$Ftot <- df$Fc+df$Fro
+
+  extreme_p <- c()
+  for(i in 1:n_tests){
+    extreme_p[i] <- min(df$exp_p[i], 1-df$exp_p[i])
+  }
+  df$extreme_p <- extreme_p
+
+  testability <- c()
+  for(i in 1:n_tests) {
+    testability[i] <- binom.test(df$Ftot[i], df$Ftot[i], df$extreme_p[i], alternative ="two.sided")$p.value
+  }
+  df$testability <- testability
+
+  # Binomial (or Chi square) Test Significance
+
+  Significance <- c()
+  for(i in 1:n_tests) {
+    ifelse(((df$Fc[i]+df$Fro[i])*(df$Ac[i]/(df$Ac[i]+df$Ao[i]))<=5 | (df$Fc[i]+df$Fro[i])*(df$Ao[i]/(df$Ac[i]+df$Ao[i]))<=5),
+           Significance[i] <- binom.test(df$Fc[i], df$Fc[i]+df$Fro[i], df$exp_p[i], alternative ="two.sided")$p.value,
+           Significance[i] <- chisq.test(c(df$Fc[i], df$Fro[i]), p = c(df$exp_p[i], 1-df$exp_p[i]))$p.value
+    )
+  }
+  df$Significance <- Significance
+
+  Test_type <- c()
+  for(i in 1:n_tests) {
+    ifelse(((df$Fc[i]+df$Fro[i])*(df$Ac[i]/(df$Ac[i]+df$Ao[i]))<=5 | (df$Fc[i]+df$Fro[i])*(df$Ao[i]/(df$Ac[i]+df$Ao[i]))<=5),
+           Test_type[i] <- "Binomial",
+           Test_type[i] <- "Chi-square"
+    )
+  }
+  df$Test_type <- Test_type
+
+  if(length(unique(df$Test_type))>1) message("Different tests were used for different canopy-recruit pairs. Check column Test_type")
+
+  Effect_int <- c()
+  for(i in 1:n_tests) {
+    ifelse((df$testability[i]>0.05),
+           Effect_int[i] <- "Not testable",
+           ifelse(df$Significance[i] > 0.05,
+                  Effect_int[i] <- "Neutral",
+                  ifelse((df$Fc[i]/df$Ac[i])>(df$Fro[i]/df$Ao[i]),
+                         Effect_int[i] <- "Facilitative",
+                         Effect_int[i] <- "Depressive")
+           )
+    )
+  }
+
+  df$Canopy_effect <- Effect_int
+  drops <- c("exp_p", "Ftot", "extreme_p")
+  df <- df[ , !(names(df) %in% drops)]
+  return(df)
+}
+
diff --git a/R/comm_to_RN.R b/R/comm_to_RN.R
new file mode 100644
index 0000000..aa1b89e
--- /dev/null
+++ b/R/comm_to_RN.R
@@ -0,0 +1,56 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param expand
+#' @param rm_sp_no_cover
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+#' # comm_to_RN: summarizes data collected by plot to the community level ( merging plots)
+#expand:expands non observed interctions
+#rm_sp_no_cover: "allsp" removes any species without cover data, while "onlycanopy":removes only canopy species without cover
+
+#that defines internaly the combination of those arguments usually used to work with different interaction types
+#rec: expands non observed interctions and removes any species without cover data (canopy or recruit)
+#fac: considers only observed intearctions and removes only canopy species without cover data (keeping recruit species without cover)
+#comp: expand non observed interactions and removes canopy species without cover data
+
+comm_to_RN<- function(int_data,cover_data, expand=c("yes","no"), rm_sp_no_cover=c("allsp","onlycanopy")){
+
+
+  if(expand=="yes" & rm_sp_no_cover=="allsp"){int_type ="rec"}
+  if(expand=="yes" & rm_sp_no_cover=="onlycanopy"){int_type ="comp"}
+  if(expand=="no" & rm_sp_no_cover=="onlycanopy"){int_type ="fac"}
+  if(expand=="no" & rm_sp_no_cover=="allsp"){int_type ="unused"}
+
+
+  if(int_type=="rec"){
+    df<-comm_to_RN_UNI(int_data,cover_data)
+  }
+
+  if(int_type=="fac"){
+    df<-comm_to_RN_BI(int_data,cover_data)
+  }
+
+  if(int_type=="comp"){
+    df<-comm_to_RN_UNI_COMP(int_data,cover_data)
+  }
+
+
+  if(int_type=="unused"){
+
+    df<-NULL
+    warning("The object is not provided because the combination of arguments is not commonly used in plant-plant interaction networks")
+  }
+
+  return(df)
+
+
+}
+
+
diff --git a/R/cum_values.R b/R/cum_values.R
new file mode 100644
index 0000000..5a71281
--- /dev/null
+++ b/R/cum_values.R
@@ -0,0 +1,87 @@
+#' Title
+#'
+#' @param int_data
+#' @param property
+#' @param k
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' #For study sites sampled in multiple plots, we can use plot accumulation curves to assess whether our estimates of these dimensions are stable. The function cum_values does plot accumulation #curves for parameters that can be obtained with package igraph. For example, for the number of #nodes and links, and connectance we can use, respectively, the igraph functions "vcount", "ecount" and "edge_density". The function may take long time if *k* >> 100 (for speed, we will run this #example with just 20 resamplings).
+# cum_values identical to cum_values_UNI
+#'
+cum_values <- function(int_data, property=c("vcount","ecount","edge_density"), k = 100){
+  require(ggplot2)
+
+  if (!"Plot" %in% names(int_data)) stop("ERROR: your interactions data lacks a column named Plots. This function requires data assembled in plots.")
+  nPlots <- length(unique(int_data$Plot))
+
+  if (nPlots < 10)
+    warning(
+      "WARNING: your are using the incidence approach with very few plots."
+    )
+
+  if(property=="vcount"){
+
+    part_RNs <- partial_RNs_UNI(int_data, k)
+    nSteps <- length(part_RNs)
+    borrar <- unlist(part_RNs, recursive = FALSE)
+    df <- data.frame(unlist(lapply(borrar, property)))
+    colnames(df) <- c("Value")
+    df$sampleSize <- sort(rep(c(1:nSteps),k))
+    plot_cumm_value <- ggplot(df, aes(x=as.factor(sampleSize), y=Value)) +
+      geom_jitter(colour="turquoise3", alpha=0.5, height = 0, width=0.1) +
+      geom_point(stat="summary", fun="mean") +
+      geom_errorbar(stat="summary", fun.data="mean_se", fun.args = list(mult = 1.96), width=0.3) +
+      labs(x="Sample Size (Num. Plots)", y="Value (mean + 95%CI)") +
+      ggtitle("Number of species")
+
+    outputs <- list("Data" = df, "Plot" = plot_cumm_value)
+    return(outputs)
+
+  }
+
+  if(property=="ecount"){
+
+    part_RNs <- partial_RNs_UNI(int_data, k)
+    nSteps <- length(part_RNs)
+    borrar <- unlist(part_RNs, recursive = FALSE)
+    df <- data.frame(unlist(lapply(borrar, property)))
+    colnames(df) <- c("Value")
+    df$sampleSize <- sort(rep(c(1:nSteps),k))
+    plot_cumm_value <- ggplot(df, aes(x=as.factor(sampleSize), y=Value)) +
+      geom_jitter(colour="turquoise3", alpha=0.5, height = 0, width=0.1) +
+      geom_point(stat="summary", fun="mean") +
+      geom_errorbar(stat="summary", fun.data="mean_se", fun.args = list(mult = 1.96), width=0.3) +
+      labs(x="Sample Size (Num. Plots)", y="Value (mean + 95%CI)") +
+      ggtitle("Number of interactions")
+
+    outputs <- list("Data" = df, "Plot" = plot_cumm_value)
+    return(outputs)
+
+  }
+
+  if(property=="edge_density"){
+
+    part_RNs <- partial_RNs_UNI(int_data, k)
+    nSteps <- length(part_RNs)
+    borrar <- unlist(part_RNs, recursive = FALSE)
+    df <- data.frame(unlist(lapply(borrar, property)))
+    colnames(df) <- c("Value")
+    df$sampleSize <- sort(rep(c(1:nSteps),k))
+    plot_cumm_value <- ggplot(df, aes(x=as.factor(sampleSize), y=Value)) +
+      geom_jitter(colour="turquoise3", alpha=0.5, height = 0, width=0.1) +
+      geom_point(stat="summary", fun="mean") +
+      geom_errorbar(stat="summary", fun.data="mean_se", fun.args = list(mult = 1.96), width=0.3) +
+      labs(x="Sample Size (Num. Plots)", y="Value (mean + 95%CI)") +
+      ggtitle("Connectance")
+
+    outputs <- list("Data" = df, "Plot" = plot_cumm_value)
+    return(outputs)
+
+  }
+
+}
+
diff --git a/R/funtopol_rec.R b/R/funtopol_rec.R
new file mode 100644
index 0000000..cdd9d5e
--- /dev/null
+++ b/R/funtopol_rec.R
@@ -0,0 +1,163 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' # antigua funtopol externa con solo la opción para redes de reclutamiento
+#'
+funtopol_rec <- function(int_data,cover_data){
+
+  int_data<-comm_to_RN(int_data,cover_data,expand="yes", rm_sp_no_cover="allsp" )
+
+  if (!"Open" %in% int_data$Canopy) stop("ERROR: your data does not contain a node named Open or it is spelled differently.")
+
+
+  int_data <- int_data[which(int_data$Fcr!=0), c("Canopy", "Recruit")]
+  g <- igraph::graph_from_data_frame(int_data, directed = TRUE)
+  g <- igraph::simplify(g, remove.multiple = TRUE, remove.loops = FALSE)
+
+  if (length(which(int_data$Canopy=="Open"))==0) {
+    warning("Open is included as a node in the network, even though no recruits (with cover data of that species) are associated with Open in this community")
+
+    g <- igraph::add_vertices(g, 1, name = "Open")
+  }
+
+  NEdges <- igraph::gsize(g)
+  NNodes <- igraph::gorder(g)
+
+  CDirected <- NEdges/(NNodes*(NNodes - 1))
+  SCCs <- igraph::components(g, mode = "strong")
+
+  if(max(SCCs$csize)>1){
+
+    numSCCs <- SCCs$no
+    numNTSCCs <- sum(SCCs$csize > 1)
+    coreSize <- max(SCCs$csize)
+    SCC_memb <- SCCs$membership
+    SCC_memb <- as.data.frame(SCC_memb)
+    SCC_subgraphs <- igraph::decompose(g, mode = "strong") # Makes a subgraph of each SCC
+    IDcore <- match(coreSize, SCCs$csize) # locates the position of the core in the list of SCCs
+    MembersCore <- igraph::V(SCC_subgraphs[[IDcore]])$name # List of the species in the core
+    IDOpen <- SCC_memb$SCC_memb[match("Open", row.names(SCC_memb))] # Locate the position of the "open" node in the list of SCCs
+    outReachFromOpen <- names(igraph::subcomponent(g, "Open", "out")[-1]) # List of nodes reachable from the "open"
+    outReachFromCore <- vector("list", coreSize) # List of nodes reachable from core nodes
+    for (i in 1:coreSize) {
+      outReachFromCore[[i]] <- igraph::subcomponent(g, MembersCore[i], mode = "out")
+    }
+    a <- unlist(outReachFromCore)
+    a <- unique(names(a))
+    MembersSatellites <- setdiff(a, MembersCore)
+    MembersTransients <- setdiff(igraph::V(g)$name,c(MembersCore,MembersSatellites))
+    MembersTransients <- MembersTransients[!MembersTransients == "Open"]
+    MembersStrictTransients <- setdiff(MembersTransients, outReachFromOpen)
+    MembersDdTransients <- setdiff(MembersTransients, MembersStrictTransients)
+    numSat <- length(MembersSatellites)
+    numTransAll <- length(MembersTransients)
+    numDdTrans <- length(MembersDdTransients)
+    numStrictTrans <- length(MembersStrictTransients)
+    propCore <- coreSize/(NNodes - 1)
+    propSat <- numSat/(NNodes - 1)
+    propTrans <- numTransAll/(NNodes - 1)
+    propStrTrans <- numStrictTrans/(NNodes - 1)
+    propDdTrans <- numDdTrans/(NNodes - 1)
+    persistence <- propCore + propSat
+
+    # Function output
+
+    df <- data.frame(
+      c(NNodes,
+        NEdges,
+        CDirected,
+        numNTSCCs,
+        coreSize,
+        propCore,
+        numSat,
+        propSat,
+        numDdTrans,
+        propDdTrans,
+        numStrictTrans,
+        propStrTrans,
+        persistence)
+    )
+    colnames(df) <- c("Value")
+    rownames(df) <- c(
+      "Num. nodes",
+      "Num. edges",
+      "Connectance",
+      "Num. non-trivial SCCs",
+      "Num. core species",
+      "Prop. core species",
+      "Num. satellite species",
+      "Prop. satellite species",
+      "Num. disturbance-dependent transients",
+      "Prop. disturbance-dependent transients",
+      "Num. strict transients",
+      "Prop. strict transients",
+      "Qualitative Persistence")
+
+    classif <- list(
+      MembersSatellites,
+      MembersCore,
+      MembersStrictTransients,
+      MembersDdTransients)
+    classif <- stats::setNames(classif,
+                               c("Satellites",
+                                 "Core",
+                                 "Strict_transients",
+                                 "Disturbance_dependent_transients")
+    )
+
+    df0 <- classif
+    df_Sat <- data.frame(df0$Satellites, rep("Satellite", length(df0$Satellites)))
+    colnames(df_Sat) <- c("id", "group")
+    df_Core <- data.frame(df0$Core, rep("Core", length(df0$Core)))
+    colnames(df_Core) <- c("id", "group")
+    df_Str <- data.frame(df0$Strict_transients, rep("Strict_transients", length(df0$Strict_transients)))
+    colnames(df_Str) <- c("id", "group")
+    df_Ddtr <- data.frame(df0$Disturbance_dependent_transients, rep("Disturbance_dependent_transients", length(df0$Disturbance_dependent_transients)))
+    colnames(df_Ddtr) <- c("id", "group")
+    df0 <- rbind(df_Sat, df_Core, df_Str, df_Ddtr)
+    df0 <- df0[order(df0$id),]
+
+    outputs <- list("Descriptors" = df, "Functional_classification" = df0)
+
+  }else{
+
+    warning("This network does not have a Core, and thus the species membership to different roles can not be defined")
+
+    coreSize<-0
+
+    df <- data.frame(
+      c(NNodes,
+        NEdges,
+        CDirected,
+        coreSize))
+
+
+    colnames(df) <- c("Value")
+    rownames(df) <- c(
+      "Num. nodes",
+      "Num. edges",
+      "Connectance",
+      "Num. core species")
+
+    df0<-list(
+      Satellites = character(0),
+      Core = character(0),
+      Strict_transients = character(0),
+      Disturbance_dependent_transients = character(0)
+    )
+
+    outputs <-list("Descriptors" = df, "Functional_classification" = df0)
+
+  }
+
+  return(outputs)
+
+}
+
diff --git a/R/int_significance.R b/R/int_significance.R
new file mode 100644
index 0000000..b8dc7e7
--- /dev/null
+++ b/R/int_significance.R
@@ -0,0 +1,46 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' # int_significance() Use all functions to consider all interactions (observed and non-observed) (for Recruitment networks) or only observed interactions (for facilitation/completition networks)
+#'
+#'
+int_significance <- function(int_data, cover_data, int_type=c("rec", "fac","comp")){
+
+  if (!"Open" %in% int_data$Canopy) stop("ERROR: tests cannot be conducted because your data does not contain a node named Open or it is spelled differently.")
+
+  if(int_type=="rec"){
+
+    data<-comm_to_RN_UNI(int_data,cover_data)
+    df<-int_significance_UNI(data)
+    if(length(unique(df$Test_type))>1) message("Different tests were used for different canopy-recruit pairs. Check column Test_type")
+  }
+
+  if(int_type=="fac"){
+
+    data<-comm_to_RN_BI(int_data,cover_data)
+    df<-int_significance_BI (data)
+    df<-df[df$Effect_int=="Enhancing",]
+    if(length(unique(df$Test_type))>1) message("Different tests were used for different canopy-recruit pairs. Check column Test_type")
+  }
+
+  if(int_type=="comp"){
+
+    data<-comm_to_RN_UNI_COMP(int_data,cover_data)
+    df<-int_significance_BI_COMP(data)
+    df<-df[df$Effect_int=="Depressing",]
+    if(length(unique(df$Test_type))>1) message("Different tests were used for different canopy-recruit pairs. Check column Test_type")
+
+  }
+
+  return(df)
+}
+
+
diff --git a/R/node_degrees.R b/R/node_degrees.R
new file mode 100644
index 0000000..e86c2a7
--- /dev/null
+++ b/R/node_degrees.R
@@ -0,0 +1,42 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' #In recruitment networks, node degrees have particular interpretations. The in- and out-degrees of #a node in a recruitment network inform, respectively, about the width of their **recruitment #niche** (number of canopy species that allow their recruitment) and of the **canopy service** they #provide (number of species that recruit in their vicinity). When the degrees are weighted by the #frequency of recruitment, then they can be interpreted, respectively, as the **abundance in the #recruit bank** (number of recruits of a species in the study site) and the **contribution of the #canopy species to the multispecific sapling bank** (number of recruits of any species associated #with the canopy species). To take into account the dominance of certain interactions in the #recruit bank, the weighted degrees can be transformed to the effective number of partners, which #can be interpreted, respectively, as the **effective recruitment niche width** and **effective #canopy service**.
+
+#desde la matriz de interacciones, de reclutamiento, facilitacion o competencia
+#para cada fila (recruit species) y columna (canopy) se cuenta la frecuencia (i.e. numero de reclutas) (service)
+#luego se hacen las matrices binarias y se cuenta el numero de especies (width)
+#'
+node_degrees <- function(int_data,cover_data, int_type=c("rec", "fac","comp")) {
+
+  if(int_type=="rec"){
+
+    data<-comm_to_RN_UNI(int_data, cover_data)
+    node_deg<-node_degrees_UNI(data)
+
+  }
+
+
+  if(int_type=="fac"){
+
+    node_deg<-suppressWarnings(node_degrees_BI(int_data,cover_data))
+
+  }
+
+  if(int_type=="comp"){
+
+    node_deg<-suppressWarnings(node_degrees_BI_COMP(int_data,cover_data))
+
+  }
+
+  return(node_deg)
+}
+
diff --git a/R/node_topol.R b/R/node_topol.R
new file mode 100644
index 0000000..01d40af
--- /dev/null
+++ b/R/node_topol.R
@@ -0,0 +1,84 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' # node_topol_UNI()node_topol(Ventisquero_com,Ventisquero_cov, int_type="rec")
+#'
+node_topol <- function(int_data,cover_data, int_type=c("rec","fac","comp")){
+
+  if(int_type=="rec"){
+
+    int_data<-comm_to_RN_UNI(int_data,cover_data)
+
+    RN_igraph <- igraph::graph_from_adjacency_matrix(t(RN_to_matrix_UNI(int_data, weight = "Pcr")), mode = "directed")
+    eigen_cent <- igraph::eigen_centrality(RN_igraph, directed=TRUE, scale=FALSE, options = list(which="LR"))$vector
+    out_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="out", mindist=1)
+    in_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="in", mindist=1)
+    df <- data.frame(eigen_cent, out_neigh,in_neigh)
+    df[, 1] <- round(df[, 1], digits = 4)
+    colnames(df) <- c("Eigenvector centrality", "Extended canopy service", "Extended recruitment niche")
+    return(df)
+
+  }
+
+  if(int_type=="fac"){
+
+    my_mat<-RN_to_matrix (int_data,cover_data, int_type="fac", weight="Pcr")
+
+    #complete with 0 to make a square matrix
+
+    all_species <- union(rownames(my_mat), colnames(my_mat))
+    new_mat <- matrix(0, nrow = length(all_species), ncol = length(all_species),
+                      dimnames = list(all_species, all_species))
+    new_mat[rownames(my_mat), colnames(my_mat)] <- my_mat
+    my_mat<-new_mat
+
+    RN_igraph <- igraph::graph_from_adjacency_matrix(t(my_mat), mode = "directed")
+    RN_igraphN <- igraph::graph_from_adjacency_matrix(my_mat, mode = "directed")
+
+    eigen_cent <- igraph::eigen_centrality(RN_igraphN, directed=TRUE, scale=FALSE, options = list(which="LR"))$vector
+    out_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="out", mindist=1)
+    in_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="in", mindist=1)
+    df <- data.frame(eigen_cent, out_neigh,in_neigh)
+    df[, 1] <- round(df[, 1], digits = 4)
+    colnames(df) <- c("Eigenvector nurse centrality", "Extended nurse service", "Extended facilitated niche")
+    return(df)
+
+
+  }
+
+  if(int_type=="comp"){
+
+    my_mat<-RN_to_matrix (int_data,cover_data, int_type="comp", weight="RII")
+    my_mat<-ifelse(my_mat<0,1,0)
+
+    #complete with 0 to make a square matrix
+    all_species <- union(rownames(my_mat), colnames(my_mat))
+    new_mat <- matrix(0, nrow = length(all_species), ncol = length(all_species),
+                      dimnames = list(all_species, all_species))
+    new_mat[rownames(my_mat), colnames(my_mat)] <- my_mat
+    my_mat<-new_mat
+
+    RN_igraph <- igraph::graph_from_adjacency_matrix(t(my_mat), mode = "directed")
+    RN_igraphN <- igraph::graph_from_adjacency_matrix(my_mat, mode = "directed")
+
+    eigen_cent <- igraph::eigen_centrality(RN_igraphN, directed=TRUE, scale=FALSE, options = list(which="LR"))$vector
+    out_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="out", mindist=1)
+    in_neigh <- igraph::neighborhood_size(RN_igraph, order=gorder(RN_igraph), mode="in", mindist=1)
+    df <- data.frame(eigen_cent, out_neigh,in_neigh)
+    df[, 1] <- round(df[, 1], digits = 4)
+    colnames(df) <- c("Eigenvector canopy centrality", "Extended canopy depression effect", "Extended recruitment depression ")
+    return(df)
+
+
+
+  }
+}
+
diff --git a/R/recruitment_niche_test.R b/R/recruitment_niche_test.R
new file mode 100644
index 0000000..b31a8b6
--- /dev/null
+++ b/R/recruitment_niche_test.R
@@ -0,0 +1,89 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' #como recruitment_niche_test_UNI pero sin quitar especies de reclutas
+#que no tienen cover ( solo las canopy que no tienen cover)
+#sustituimos pre_associndex_UNISITE_UNI por
+#associndex(int_data,cover_data, expand="yes", rm_sp_no_cover="onlycanopy")
+#calcula para cada especie de recluta, si se establece con una mayor (o menor)
+# densidad de reclutas (bajo cualquier especie de canopy) que la esperada en funcion del area que ocupa las canopies en conjunto y el open.
+#The input data are two files: the field data containing interactions and species cover
+#'
+recruitment_niche_test <- function(int_data,cover_data){
+
+  df <- associndex(int_data,cover_data, expand="yes", rm_sp_no_cover="onlycanopy")
+  sp_Fr <- stats::aggregate(Fcr ~ Recruit, data = df, FUN = sum)
+  sp_Av <- stats::aggregate(Ac ~ Recruit, data = df, FUN = sum)
+  sp_Fro <- stats::aggregate(Fro ~ Recruit, data = df, FUN = max)
+  sp_Ao <- stats::aggregate(Ao ~ Recruit, data = df, FUN = max)
+  n_tests <- dim(sp_Fr)[1]
+  df <- data.frame(c(sp_Fr, sp_Av, sp_Fro, sp_Ao))
+  myvars <- names(df) %in% c("Recruit.1", "Recruit.2", "Recruit.3")
+  df <- df[!myvars]
+  colnames(df) <- c("Recruit", "Fr", "Av", "Fro", "Ao")
+  df$exp_p <- df$Av/(df$Av+df$Ao) # Expected probability of success (i.e. of recruiting under canopy)
+
+  # Testability through Binomial test
+
+  df$Ftot <- df$Fr+df$Fro
+
+  extreme_p <- c()
+  for(i in 1:n_tests){
+    extreme_p[i] <- min(df$exp_p[i], 1-df$exp_p[i])
+  }
+  df$extreme_p <- extreme_p
+
+  testability <- c()
+  for(i in 1:n_tests) {
+    testability[i] <- binom.test(df$Ftot[i], df$Ftot[i], df$extreme_p[i], alternative ="two.sided")$p.value
+  }
+  df$testability <- testability
+
+  # Binomial (or Chi square) Test Significance
+
+  Significance <- c()
+  for(i in 1:n_tests) {
+    ifelse(((df$Fr[i]+df$Fro[i])*(df$Av[i]/(df$Av[i]+df$Ao[i]))<=5 | (df$Fr[i]+df$Fro[i])*(df$Ao[i]/(df$Av[i]+df$Ao[i]))<=5),
+           Significance[i] <- binom.test(df$Fr[i], df$Fr[i]+df$Fro[i], df$exp_p[i], alternative ="two.sided")$p.value,
+           Significance[i] <- chisq.test(c(df$Fr[i], df$Fro[i]), p = c(df$exp_p[i], 1-df$exp_p[i]))$p.value
+    )
+  }
+  df$Significance <- Significance
+
+  Test_type <- c()
+  for(i in 1:n_tests) {
+    ifelse(((df$Fr[i]+df$Fro[i])*(df$Av[i]/(df$Av[i]+df$Ao[i]))<=5 | (df$Fr[i]+df$Fro[i])*(df$Ao[i]/(df$Av[i]+df$Ao[i]))<=5),
+           Test_type[i] <- "Binomial",
+           Test_type[i] <- "Chi-square"
+    )
+  }
+  df$Test_type <- Test_type
+
+  if(length(unique(df$Test_type))>1) warning("Different tests were used for different canopy-recruit pairs. Check column Test_type")
+
+  Effect_int <- c()
+  for(i in 1:n_tests) {
+    ifelse((df$testability[i]>0.05),
+           Effect_int[i] <- "Not testable",
+           ifelse(df$Significance[i] > 0.05,
+                  Effect_int[i] <- "Neutral",
+                  ifelse((df$Fr[i]/df$Av[i])>(df$Fro[i]/df$Ao[i]),
+                         Effect_int[i] <- "Facilitated",
+                         Effect_int[i] <- "Depressed")
+           )
+    )
+  }
+
+  df$Veg_effect <- Effect_int
+  drops <- c("exp_p", "Ftot", "extreme_p")
+  df <- df[ , !(names(df) %in% drops)]
+  return(df)
+}
+
diff --git a/R/remove_no_cover.R b/R/remove_no_cover.R
new file mode 100644
index 0000000..38ef1e0
--- /dev/null
+++ b/R/remove_no_cover.R
@@ -0,0 +1,33 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param rm_sp_no_cover
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+# remove_no_cover : remove the species from the interaction database that do not have cover data
+#rm_sp_no_cover:
+#"allsp" removes any species without cover data,
+#"onlycanopy":removes only canopy species without cover
+
+#for rec interaction type remove any species without cover data from
+#the interactions database, and for comp or fac remove only canopy species without cover data
+#keeping recruit speceis without cover data
+#'
+remove_no_cover<- function(int_data,cover_data, rm_sp_no_cover=c("allsp","onlycanopy")){
+
+  if(rm_sp_no_cover=="allsp"){
+    df<-remove_no_cover_UNI(int_data,cover_data)
+  }
+
+  if(rm_sp_no_cover=="onlycanopy"){
+    df<-remove_no_cover_BI(int_data,cover_data)
+  }
+
+  return(df)
+}
diff --git a/R/topol_depre.R b/R/topol_depre.R
new file mode 100644
index 0000000..fa15b1d
--- /dev/null
+++ b/R/topol_depre.R
@@ -0,0 +1,336 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param direction
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+topol_depre <- function(int_data,cover_data, direction=c("in","out")){
+
+  direction <- match.arg(direction)
+
+  M<-RN_to_matrix(int_data, cover_data, int_type="comp",weight="Pcr")
+
+  #make a square matrix to build a unipartite directed graph to assess reciprocal facilitation
+
+  species <- union(rownames(M), colnames(M))
+
+  A <- matrix(0,  nrow = length(species),ncol = length(species),dimnames = list(species, species))
+  A[rownames(M), colnames(M)] <- M
+
+  # traspose the matrix to set that species in rows (from) facilitate species in columns (to)
+
+  M<-t(A)
+
+  #generate a graph (A facilitates B)
+
+  g<-graph_from_adjacency_matrix(M, mode="directed", diag=FALSE)
+
+  ###generate a list to save the objects
+
+  indirect_comp <- list(
+    loops = list(
+      summary = list(),
+      nodes   = list()
+    ),
+    simple = list(
+      summary = list(),
+      nodes   = list()
+    )
+  )
+
+  ######
+  #Loops of reciprocal facilitation
+  #######
+
+  #nodes
+  scc <- components(g, mode = "strong")
+  scc_groups <- split(V(g)$name, scc$membership)
+
+  scc_groups<-scc_groups[sapply(scc_groups, length) > 1]
+
+  #summary
+  table_scc <- data.frame(
+    scc_id = seq_along(scc_groups),
+    n_nodos = sapply(scc_groups, length)
+  )
+
+
+  ########
+  #functions to detect simple paths and filter subpaths
+  ########
+
+  #function to detect subpaths
+
+  es_subpath <- function(p, q) {
+    lp <- length(p)
+    lq <- length(q)
+
+    if (lp >= lq) return(FALSE)
+
+    for (i in 1:(lq - lp + 1)) {
+      if (all(q[i:(i + lp - 1)] == p)) {
+        return(TRUE)
+      }
+    }
+    FALSE
+  }
+
+  #function to filter subpaths
+
+  filter_paths <- function(paths) {
+
+    lens <- sapply(paths, length)
+    ord <- order(lens, decreasing = TRUE)
+    paths <- paths[ord]
+
+    keep <- rep(TRUE, length(paths))
+
+    for (i in seq_along(paths)) {
+
+      if (!keep[i]) next
+
+      # compare with paths already aceptaded (& longer)
+      if (i > 1) {
+        for (j in which(keep)[which(keep) < i]) {
+          if (es_subpath(paths[[i]], paths[[j]])) {
+            keep[i] <- FALSE
+            break
+          }
+        }
+      }
+    }
+
+    paths[keep]
+  }
+
+
+
+
+  #############
+  ##Function to obtain simple paths (linear)
+  #############
+
+  #function to get maximal paths
+
+
+
+  #the threshodl indicates the percentage of overlap of two paths to me considered distinct
+  maximal_paths_distinct_pruned <- function(g, start, overlap_threshold = 0.75, dir=c("in","out")) {
+
+    best_paths <- list()  # paths filtered by overlap
+
+    #Internal function applying Depth-First Search (DFS) algorithm and prunning
+    #How does it work: DFS starts at a root node.
+    #It explores an unvisited neighbor and goes “deep” until it can’t go any further.
+    #Then it backtracks to explore other neighbors.
+    #This process repeats until all reachable nodes have been visited.
+
+    dfs_extend <- function(current_path) {
+
+      last_node <- tail(current_path, 1)
+
+      if(dir=="out"){
+        neighbors_exp <- setdiff(neighbors(g, last_node, mode = "out")$name, current_path)
+
+      }
+
+      if(dir=="in"){
+        neighbors_exp <- setdiff(neighbors(g, last_node, mode = "in")$name, current_path)
+
+      }
+      # Basic case: it can not be extended
+      if (length(neighbors_exp) == 0) {
+
+        keep <- TRUE
+        remove_idx <- c()
+
+        for (i in seq_along(best_paths)) {
+          existing <- best_paths[[i]]
+
+          # calculate oeverlap between founded paths
+          n_shared <- length(intersect(current_path, existing))
+          overlap <- n_shared / min(length(current_path), length(existing))
+
+          if (overlap >= overlap_threshold) {
+            # path too similar to an existing one → discard
+            keep <- FALSE
+            break
+          }
+
+          # remove existing paths that are included in the new one
+          if (all(existing %in% current_path)) {
+            remove_idx <- c(remove_idx, i)
+          }
+        }
+
+        if (keep) {
+          if (length(remove_idx) > 0) best_paths[remove_idx] <<- NULL
+          best_paths[[length(best_paths) + 1]] <<- current_path
+        }
+
+        return()
+      }
+
+      # early prunning: if the potential extension is going to produce a lot of overlap
+      for (n in neighbors_exp) {
+        potential_path <- c(current_path, n)
+
+        prune <- FALSE
+        for (existing in best_paths) {
+          n_shared <- length(intersect(potential_path, existing))
+          overlap <- n_shared / min(length(potential_path), length(existing))
+          if (overlap >= overlap_threshold) {
+            prune <- TRUE
+            break
+          }
+        }
+
+        if (!prune) {
+          dfs_extend(potential_path)
+        }
+      }
+    }
+
+    dfs_extend(start)
+    return(best_paths)
+  }
+
+
+  if(direction=="out"){
+
+    #############
+    ##Obtain simple paths (linear) FROM all nodes
+    #############
+
+
+    #nodes
+
+    # list to save the linea paths for each node
+
+    paths_per_node_out <- vector("list", vcount(g))
+
+    for (i in seq_along(V(g))) {
+
+      nodo <- V(g)[i]
+
+      # subgrafo reachable from this node
+      reachable <- subcomponent(g, nodo, mode = "out")
+      g_sub <- induced_subgraph(g, reachable)
+
+      # eliminate edges that come back to the initial node
+      edges_to_delete <- E(g_sub)[.to(V(g_sub)[name == V(g)$name[i]])]
+      g_sub <- delete_edges(g_sub, edges_to_delete)
+
+      # calculate máximum paths
+      paths <- maximal_paths_distinct_pruned(g_sub, V(g_sub)$name[V(g_sub)$name == V(g)$name[i]], dir="out")
+
+      # save
+      paths_per_node_out[[i]] <- paths
+    }
+
+    names(paths_per_node_out) <- V(g)$name
+
+
+    #summary
+
+    table_paths_out <- lapply(names(paths_per_node_out), function(nodo) {
+
+      paths <- paths_per_node_out[[nodo]]
+
+      # number of paths from this node
+      n_paths <- length(paths)
+
+      # number of nodes in each path
+      path_lengths <- sapply(paths, length)
+
+      # generate a data.frame per node
+      data.frame(
+        nodo = nodo,
+        path_index = seq_len(n_paths),
+        n_nodes_in_path = path_lengths,
+        stringsAsFactors = FALSE
+      )
+    })
+
+    # combine in a single data.frame
+    table_paths_out <- do.call(rbind, table_paths_out)
+
+    #save objects
+
+    indirect_comp$loops$summary<-table_scc
+    indirect_comp$loops$nodes<-scc_groups
+    indirect_comp$simple$summary<-table_paths_out
+    indirect_comp$simple$nodes<-paths_per_node_out
+
+  }
+
+  #############
+  ##Obtain simple paths (linear) TO all nodes
+  #############
+
+  if(direction=="in"){
+
+    #nodes
+
+    paths_per_node_in <- vector("list", vcount(g))
+
+    for (i in seq_along(V(g))) {
+
+      nodo <- V(g)[i]
+
+      # subgrafo reaching this node
+      reachable <- subcomponent(g, nodo, mode = "in")
+      g_sub <- induced_subgraph(g, reachable)
+
+
+      # calculate máximum paths distinct between them
+      paths <- maximal_paths_distinct_pruned(g_sub, V(g_sub)$name[V(g_sub)$name == V(g)$name[i]], dir="in")
+
+      # save
+      paths_per_node_in[[i]] <- paths
+    }
+
+    names(paths_per_node_in) <- V(g)$name
+
+
+
+    #summary
+
+    table_paths_in <- lapply(names(paths_per_node_in), function(nodo) {
+
+      paths <- paths_per_node_in[[nodo]]
+
+      # number of paths reaching this node
+      n_paths <- length(paths)
+
+      # número de nodos en cada path
+      path_lengths <- sapply(paths, length)
+
+      # generar un data.frame por nodo
+      data.frame(
+        nodo = nodo,
+        path_index = seq_len(n_paths),
+        n_nodes_in_path = path_lengths,
+        stringsAsFactors = FALSE
+      )
+    })
+
+    # combine in a single data.frame
+    table_paths_in <- do.call(rbind, table_paths_in)
+
+    #save the objects
+    indirect_comp$loops$summary<-table_scc
+    indirect_comp$loops$nodes<-scc_groups
+    indirect_comp$simple$summary<-table_paths_in
+    indirect_comp$simple$nodes<-paths_per_node_in
+
+  }
+
+  return(indirect_comp)
+}
diff --git a/R/topol_fac.R b/R/topol_fac.R
new file mode 100644
index 0000000..3c3b06e
--- /dev/null
+++ b/R/topol_fac.R
@@ -0,0 +1,338 @@
+
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param direction
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+topol_fac <- function(int_data,cover_data, direction=c("in","out")){
+
+  direction <- match.arg(direction)
+
+  M<-RN_to_matrix(int_data, cover_data, int_type="fac",weight="Pcr")
+
+  #make a square matrix to build a unipartite directed graph to assess reciprocal facilitation
+
+  species <- union(rownames(M), colnames(M))
+
+  A <- matrix(0,  nrow = length(species),ncol = length(species),dimnames = list(species, species))
+  A[rownames(M), colnames(M)] <- M
+
+  # traspose the matrix to set that species in rows (from) facilitate species in columns (to)
+
+  M<-t(A)
+
+  #generate a graph (A facilitates B)
+
+  g<-graph_from_adjacency_matrix(M, mode="directed", diag=FALSE)
+
+  ###generate a list to save the objects
+
+  indirect_fac <- list(
+    loops = list(
+      summary = list(),
+      nodes   = list()
+    ),
+    simple = list(
+      summary = list(),
+      nodes   = list()
+    )
+  )
+
+  ######
+  #Loops of reciprocal facilitation
+  #######
+
+  #nodes
+  scc <- components(g, mode = "strong")
+  scc_groups <- split(V(g)$name, scc$membership)
+
+  scc_groups<-scc_groups[sapply(scc_groups, length) > 1]
+
+  #summary
+  table_scc <- data.frame(
+    scc_id = seq_along(scc_groups),
+    n_nodos = sapply(scc_groups, length)
+  )
+
+
+  ########
+  #functions to detect simple paths and filter subpaths
+  ########
+
+  #function to detect subpaths
+
+  es_subpath <- function(p, q) {
+    lp <- length(p)
+    lq <- length(q)
+
+    if (lp >= lq) return(FALSE)
+
+    for (i in 1:(lq - lp + 1)) {
+      if (all(q[i:(i + lp - 1)] == p)) {
+        return(TRUE)
+      }
+    }
+    FALSE
+  }
+
+  #function to filter subpaths
+
+  filter_paths <- function(paths) {
+
+    lens <- sapply(paths, length)
+    ord <- order(lens, decreasing = TRUE)
+    paths <- paths[ord]
+
+    keep <- rep(TRUE, length(paths))
+
+    for (i in seq_along(paths)) {
+
+      if (!keep[i]) next
+
+      # compare with paths already aceptaded (& longer)
+      if (i > 1) {
+        for (j in which(keep)[which(keep) < i]) {
+          if (es_subpath(paths[[i]], paths[[j]])) {
+            keep[i] <- FALSE
+            break
+          }
+        }
+      }
+    }
+
+    paths[keep]
+  }
+
+
+
+
+  #############
+  ##Funcions to obtain simple paths (linear)
+  #############
+
+
+  #function to get maximal paths
+
+
+
+  #the threshodl indicates the percentage of overlap of two paths to me considered distinct
+
+  maximal_paths_distinct_pruned <- function(g, start, overlap_threshold = 0.75, dir=c("in","out")) {
+
+    best_paths <- list()  # paths filtered by overlap
+
+    #Internal function applying Depth-First Search (DFS) algorithm and prunning
+    #How does it work: DFS starts at a root node.
+    #It explores an unvisited neighbor and goes “deep” until it can’t go any further.
+    #Then it backtracks to explore other neighbors.
+    #This process repeats until all reachable nodes have been visited.
+
+    dfs_extend <- function(current_path) {
+
+      last_node <- tail(current_path, 1)
+
+      if(dir=="out"){
+        neighbors_exp <- setdiff(neighbors(g, last_node, mode = "out")$name, current_path)
+
+      }
+
+      if(dir=="in"){
+        neighbors_exp <- setdiff(neighbors(g, last_node, mode = "in")$name, current_path)
+
+      }
+
+
+      # Basic case: it can not be extended
+      if (length(neighbors_exp) == 0) {
+
+        keep <- TRUE
+        remove_idx <- c()
+
+        for (i in seq_along(best_paths)) {
+          existing <- best_paths[[i]]
+
+          # calculate oeverlap between founded paths
+          n_shared <- length(intersect(current_path, existing))
+          overlap <- n_shared / min(length(current_path), length(existing))
+
+          if (overlap >= overlap_threshold) {
+            # path too similar to an existing one → discard
+            keep <- FALSE
+            break
+          }
+
+          # remove existing paths that are included in the new one
+          if (all(existing %in% current_path)) {
+            remove_idx <- c(remove_idx, i)
+          }
+        }
+
+        if (keep) {
+          if (length(remove_idx) > 0) best_paths[remove_idx] <<- NULL
+          best_paths[[length(best_paths) + 1]] <<- current_path
+        }
+
+        return()
+      }
+
+      # early prunning: if the potential extension is going to produce a lot of overlap
+      for (n in neighbors_exp) {
+        potential_path <- c(current_path, n)
+
+        prune <- FALSE
+        for (existing in best_paths) {
+          n_shared <- length(intersect(potential_path, existing))
+          overlap <- n_shared / min(length(potential_path), length(existing))
+          if (overlap >= overlap_threshold) {
+            prune <- TRUE
+            break
+          }
+        }
+
+        if (!prune) {
+          dfs_extend(potential_path)
+        }
+      }
+    }
+
+    dfs_extend(start)
+    return(best_paths)
+  }
+
+  if(direction=="out"){
+
+    #############
+    ##Obtain simple paths (linear) FROM all nodes
+    #############
+    #nodes
+
+    # list to save the linea paths for each node
+
+    paths_per_node_out <- vector("list", vcount(g))
+
+    for (i in seq_along(V(g))) {
+
+      nodo <- V(g)[i]
+
+      # subgrafo reachable from this node
+      reachable <- subcomponent(g, nodo, mode = "out")
+      g_sub <- induced_subgraph(g, reachable)
+
+      # eliminate edges that come back to the initial node
+      edges_to_delete <- E(g_sub)[.to(V(g_sub)[name == V(g)$name[i]])]
+      g_sub <- delete_edges(g_sub, edges_to_delete)
+
+      # calculate máximum paths
+      paths <- maximal_paths_distinct_pruned(g_sub, V(g_sub)$name[V(g_sub)$name == V(g)$name[i]],dir="out")
+
+      # save
+      paths_per_node_out[[i]] <- paths
+    }
+
+    names(paths_per_node_out) <- V(g)$name
+
+
+    #summary
+
+    table_paths_out <- lapply(names(paths_per_node_out), function(nodo) {
+
+      paths <- paths_per_node_out[[nodo]]
+
+      # number of paths from this node
+      n_paths <- length(paths)
+
+      # number of nodes in each path
+      path_lengths <- sapply(paths, length)
+
+      # generate a data.frame per node
+      data.frame(
+        nodo = nodo,
+        path_index = seq_len(n_paths),
+        n_nodes_in_path = path_lengths,
+        stringsAsFactors = FALSE
+      )
+    })
+
+    # combine in a single data.frame
+    table_paths_out <- do.call(rbind, table_paths_out)
+
+    #save objects
+
+    indirect_fac$loops$summary<-table_scc
+    indirect_fac$loops$nodes<-scc_groups
+    indirect_fac$simple$summary<-table_paths_out
+    indirect_fac$simple$nodes<-paths_per_node_out
+
+  }
+
+  #############
+  ##Obtain simple paths (linear) TO all nodes
+  #############
+
+  if(direction=="in"){
+
+    #nodes
+
+    paths_per_node_in <- vector("list", vcount(g))
+
+    for (i in seq_along(V(g))) {
+
+      nodo <- V(g)[i]
+
+      # subgrafo reaching this node
+      reachable <- subcomponent(g, nodo, mode = "in")
+      g_sub <- induced_subgraph(g, reachable)
+
+
+      # calculate máximum paths distinct between them
+      paths <- maximal_paths_distinct_pruned(g_sub, V(g_sub)$name[V(g_sub)$name == V(g)$name[i]], dir="in")
+
+      # save
+      paths_per_node_in[[i]] <- paths
+    }
+
+    names(paths_per_node_in) <- V(g)$name
+
+
+
+    #summary
+
+    table_paths_in <- lapply(names(paths_per_node_in), function(nodo) {
+
+      paths <- paths_per_node_in[[nodo]]
+
+      # number of paths reaching this node
+      n_paths <- length(paths)
+
+      # número de nodos en cada path
+      path_lengths <- sapply(paths, length)
+
+      # generar un data.frame por nodo
+      data.frame(
+        nodo = nodo,
+        path_index = seq_len(n_paths),
+        n_nodes_in_path = path_lengths,
+        stringsAsFactors = FALSE
+      )
+    })
+
+    # combine in a single data.frame
+    table_paths_in <- do.call(rbind, table_paths_in)
+
+    #save the objects
+    indirect_fac$loops$summary<-table_scc
+    indirect_fac$loops$nodes<-scc_groups
+    indirect_fac$simple$summary<-table_paths_in
+    indirect_fac$simple$nodes<-paths_per_node_in
+
+  }
+
+  return(indirect_fac)
+}
diff --git a/R/visu_funtopol_rec.R b/R/visu_funtopol_rec.R
new file mode 100644
index 0000000..0e1e316
--- /dev/null
+++ b/R/visu_funtopol_rec.R
@@ -0,0 +1,75 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' # antigua visu_funtopol externa con solo la opción para redes de reclutamiento
+#'necesita el paquete visNetwork(segun Rstudio)
+#'
+#'
+visu_funtopol_rec <- function(int_data,cover_data){
+
+  int_data0<-int_data
+  cover_data0<-cover_data
+
+  int_data<-comm_to_RN(int_data,cover_data,expand="yes", rm_sp_no_cover="allsp" )
+
+  if (!"Open" %in% int_data$Canopy) stop("ERROR: your data does not contain a node named Open or it is spelled differently.")
+
+
+  nodes_list <- funtopol_UNI(int_data)$Functional_classification
+
+  int_data_plot<-int_significance(int_data0,cover_data0, int_type=c("rec"))
+  int_data_plot <- int_data_plot[, c("Canopy", "Recruit", setdiff(names(int_data_plot), c("Canopy","Recruit")))]
+
+  g <- igraph::graph_from_data_frame(int_data_plot, directed = TRUE)
+  g <- igraph::simplify(g, remove.multiple = TRUE, remove.loops = FALSE)
+  SCCs <- igraph::components(g, mode = "strong")
+
+
+  core<-funtopol_UNI(int_data)$Descriptors[which(rownames(funtopol_UNI(int_data)$Descriptors)=="Num. core species"),]
+  if(core>0){
+
+    open_df <- c("Open", "Open")
+    nodes_list <- rbind(nodes_list, open_df)
+    nodes_list$label <- nodes_list$id
+    int_data <- int_data[which(int_data$Fcr!=0), c("Canopy", "Recruit")]
+    g <- igraph::graph_from_data_frame(int_data, directed = TRUE)
+    g <- igraph::simplify(g, remove.multiple = TRUE, remove.loops = FALSE)
+    g<- igraph::as_data_frame(g, what = "both")
+    edges_list <- g$edges
+
+    # nodes data.frame for legend
+    lnodes <- data.frame(label = c("Open", "Core", "Satellite", "Strict transient", "Disturbance-dependent transient"),
+                         shape = c( "dot"), color = c("#F0E442", "#009E73", "#0072B2", "#D55E00", "#CC79A7"),
+                         title = "Functional types", id = 1:5)
+
+    # Network visualization and export to html
+    require(visNetwork)
+    network <- visNetwork(nodes_list, edges_list) %>%
+      visNetwork::visIgraphLayout(layout = "layout_with_fr") %>%
+      visEdges(arrows ="to") %>%
+      visGroups(groupname = "Open", color = "#F0E442") %>%
+      visGroups(groupname = "Core", color = "#009E73") %>%
+      visGroups(groupname = "Satellite", color = "#0072B2") %>%
+      visGroups(groupname = "Strict_transients", color = "#D55E00") %>%
+      visGroups(groupname = "Disturbance_dependent_transients", color = "#CC79A7") %>%
+      visOptions(nodesIdSelection = TRUE) %>%
+      visNetwork::visLegend(addNodes = lnodes, useGroups = FALSE)
+
+
+    htmlwidgets::saveWidget(network, file = "network_recruitment.html", selfcontained = FALSE)# Save the html version of the network
+
+  }else{
+
+    stop("ERROR: This network does not have a Core, and thus the functional topology can not be visualized")}
+
+
+  return(network)
+}
+
diff --git a/R/visu_net.R b/R/visu_net.R
new file mode 100644
index 0000000..9bb991b
--- /dev/null
+++ b/R/visu_net.R
@@ -0,0 +1,180 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param int_type
+#' @param weight
+#' @param mode
+#' @param scale_w
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#' #Nueva funcion para visualizar la red de interacciones
+#'
+visu_net<-function(int_data,cover_data,int_type=c("rec","fac","comp"), weight = c("Pcr","Fcr","Dcr","Dro","Ns", "NintC", "NintA", "RII"), mode= c("uni","bi"), scale_w=1) {
+
+
+
+  int_type <- match.arg(int_type)
+  weight <- match.arg(weight)
+  mode <- match.arg(mode)
+
+  if(int_type=="rec" && mode=="bi"){
+    warning("Here recruitment networks are considered as replacement networks sensu Alcantara et al. 2019. JVS, 30:1239-1249, and thus unipartite by definition. Therefore, the combination of int_type=rec,and mode=uni is not provided")
+  }
+
+  if(int_type=="rec" && mode=="uni"){
+
+    mat <- suppressWarnings(t(RN_to_matrix(int_data, cover_data, int_type = "rec", weight = weight)))
+    edge_list <- as.data.frame(as.table(mat))
+    colnames(edge_list) <- c("from", "to", "weight")
+    edge_list <- subset(edge_list, weight > 0)
+    RN_igraph <- graph_from_data_frame(edge_list, directed = TRUE)
+
+    if (weight %in% c("Ns", "NintC", "NintA", "RII")) {
+      stop("the index specified in the weight argument uses open as a reference, meanwhile recruitment networks  considered it as a node within the network. This creates an inconsistency as open cannot simultaneously function as a node in the network and as a baseline for weighting interactions.")
+    }
+
+
+    scale<-scale_w # a factor used to adjust the magnitude of the weight ( i.e. width of the arrows)
+    # We transpose the adjacency matrix so that arrows point from canopy to recruit,
+    #this represents which species will replace a the space ocupied by the canopy in the future.
+    plot(RN_igraph,
+         edge.arrow.size=.3,
+         edge.width = E(RN_igraph)$weight*scale,
+         vertex.color="chartreuse",
+         vertex.size=8,
+         vertex.frame.color="darkolivegreen",
+         vertex.label.color="black",
+         vertex.label.cex=0.8,
+         vertex.label.dist=2,
+         vertex.label.font = 3,
+         edge.curved=0.2,
+         #layout=layout_with_kk(RN_igraph),
+         layout=layout_in_circle(RN_igraph),
+         frame = TRUE)
+    title(main="Recruitment Network")
+    return(RN_igraph)
+
+
+
+  }
+
+
+  if(int_type=="fac" && mode=="uni"){
+
+
+    mat <- RN_to_matrix(int_data, cover_data, int_type = "fac", weight = weight)
+    edge_list <- as.data.frame(as.table(mat))
+    colnames(edge_list) <- c("from", "to", "weight")
+    edge_list <- subset(edge_list, weight > 0)
+
+    RN_igraph <- graph_from_data_frame(edge_list, directed = TRUE)
+
+    scale<-scale_w # a factor used to adjust the magnitude of the weight ( i.e. width of the arrows)
+    # We transpose the adjacency matrix so that arrows point from canopy to recruit,
+    #this represents which species enhances its recuitment under another species or itself.
+    plot(RN_igraph,
+         edge.arrow.size=.3,
+         edge.width = E(RN_igraph)$weight*scale,
+         vertex.color="chartreuse",
+         vertex.size=8,
+         vertex.frame.color="darkolivegreen",
+         vertex.label.color="black",
+         vertex.label.cex=0.8,
+         vertex.label.dist=2,
+         vertex.label.font = 3,
+         edge.curved=0.2,
+         #layout=layout_with_kk(RN_igraph),
+         layout=layout_in_circle(RN_igraph),
+         frame = TRUE)
+    title(main="Unipartite Recruitment Enhancement Network")
+
+    return(RN_igraph)
+  }
+
+
+  if(int_type=="comp" && mode=="uni"){
+
+
+    mat <- RN_to_matrix(int_data, cover_data, int_type = "comp", weight = weight)
+    edge_list <- as.data.frame(as.table(mat))
+    colnames(edge_list) <- c("from", "to", "weight")
+
+    if (weight %in% c("Ns", "NintC", "NintA", "RII")) {
+      edge_list$weight<-abs(edge_list$weight)
+    }
+    edge_list <- subset(edge_list, weight > 0)
+    RN_igraph <- graph_from_data_frame(edge_list, directed = TRUE)
+
+
+    scale<-scale_w # a factor used to adjust the magnitude of the weight ( i.e. width of the arrows)
+    # We transpose the adjacency matrix so that arrows point from canopy to recruit,
+    #this represents which species enhances its recuitment under another species or itself.
+    plot(RN_igraph,
+         edge.arrow.size=.3,
+         edge.width = E(RN_igraph)$weight*scale,
+         vertex.color="chartreuse",
+         vertex.size=8,
+         vertex.frame.color="darkolivegreen",
+         vertex.label.color="black",
+         vertex.label.cex=0.8,
+         vertex.label.dist=2,
+         vertex.label.font = 3,
+         edge.curved=0.2,
+         #layout=layout_with_kk(RN_igraph),
+         layout=layout_in_circle(RN_igraph),
+         frame = TRUE)
+    title(main="Unipartite Recruitment Depression Network")
+
+    return(RN_igraph)
+
+  }
+
+
+  if(int_type=="fac" && mode=="bi"){
+
+
+    a <- RN_to_matrix(int_data, cover_data, int_type = "fac", weight = weight)
+    scale<-scale_w
+    a<-a*scale
+    sorted_a <- sortweb(a, sort.order = "dec")
+
+    # Graficar la red
+    plotweb(sorted_a,
+            srt = 90,
+            higher_italic = TRUE,
+            lower_italic = TRUE,
+            higher_color = "#E69F00",  # Canopy
+            lower_color = "#0072B2"   # Recruit
+    )
+
+
+  }
+
+
+  if(int_type=="comp" && mode=="bi"){
+
+
+    a <- RN_to_matrix(int_data, cover_data, int_type = "comp", weight = weight)
+    scale<-scale_w
+
+    if(weight%in%c("Ns", "NintC", "NintA", "RII")){a<-a*scale*(-1)}else{a<-a*scale}
+
+    sorted_a <- sortweb(a, sort.order = "dec")
+
+    # Graficar la red
+    plotweb(sorted_a,
+            srt = 90,
+            higher_italic = TRUE,
+            lower_italic = TRUE,
+            higher_color = "#E69F00",  # Canopy
+            lower_color = "#0072B2"   # Recruit
+    )
+
+  }
+}
+
diff --git a/R/visu_topol_depre.R b/R/visu_topol_depre.R
new file mode 100644
index 0000000..d81f478
--- /dev/null
+++ b/R/visu_topol_depre.R
@@ -0,0 +1,72 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param layout_fun
+#' @param vertex_size
+#' @param edge_arrow_size
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+visu_topol_depre <- function(int_data,
+                             cover_data,
+                             layout_fun = layout_with_fr,
+                             vertex_size = 20,
+                             edge_arrow_size = 0.4) {
+
+
+  # 1. Rebuild graph (same logic as inside depre_topol)
+  M <- RN_to_matrix(int_data, cover_data, int_type = "comp", weight = "Pcr")
+
+  species <- union(rownames(M), colnames(M))
+
+  A <- matrix(0,
+              nrow = length(species),
+              ncol = length(species),
+              dimnames = list(species, species))
+
+  A[rownames(M), colnames(M)] <- M
+  M <- t(A)
+
+  g <- graph_from_adjacency_matrix(M,
+                                   mode = "directed",
+                                   diag = FALSE)
+
+  # 2. Detect SCCs
+  scc <- igraph::components(g, mode = "strong")
+  scc_groups <- split(V(g)$name, scc$membership)
+
+  # Keep only real loops (>1 node)
+  scc_groups <- scc_groups[sapply(scc_groups, length) > 1]
+
+  # 3. Assign colors
+  vertex_colors <- rep("grey80", vcount(g))
+  names(vertex_colors) <- V(g)$name
+
+  if (length(scc_groups) > 0) {
+
+    palette_colors <- rainbow(length(scc_groups))
+
+    for (i in seq_along(scc_groups)) {
+      vertex_colors[scc_groups[[i]]] <- palette_colors[i]
+    }
+  }
+
+  # 4. Plot
+  plot(g,
+       layout = layout_fun(g),
+       vertex.color = vertex_colors,
+       vertex.size = vertex_size,
+       vertex.label.cex = 0.8,
+       vertex.frame.color = "black",
+       edge.arrow.size = edge_arrow_size,
+       main = paste0("Depression recruitment loops (Canopy -> Recruit)"))
+
+  # Return invisibly for further manipulation if needed
+  invisible(list(graph = g,
+                 scc = scc_groups))
+}
diff --git a/R/visu_topol_fac.R b/R/visu_topol_fac.R
new file mode 100644
index 0000000..79b60d9
--- /dev/null
+++ b/R/visu_topol_fac.R
@@ -0,0 +1,74 @@
+#' Title
+#'
+#' @param int_data
+#' @param cover_data
+#' @param layout_fun
+#' @param vertex_size
+#' @param edge_arrow_size
+#'
+#' @returns
+#' @export
+#'
+#' @examples
+#'
+#'
+visu_topol_fac <- function(int_data,
+                           cover_data,
+                           layout_fun = layout_with_fr,
+                           vertex_size = 20,
+                           edge_arrow_size = 0.4) {
+
+
+
+
+  # 1. Rebuild graph (same logic as inside depre_topol)
+  M <- RN_to_matrix(int_data, cover_data, int_type = "fac", weight = "Pcr")
+
+  species <- union(rownames(M), colnames(M))
+
+  A <- matrix(0,
+              nrow = length(species),
+              ncol = length(species),
+              dimnames = list(species, species))
+
+  A[rownames(M), colnames(M)] <- M
+  M <- t(A)
+
+  g <- graph_from_adjacency_matrix(M,
+                                   mode = "directed",
+                                   diag = FALSE)
+
+  # 2. Detect SCCs
+  scc <- igraph::components(g, mode = "strong")
+  scc_groups <- split(V(g)$name, scc$membership)
+
+  # Keep only real loops (>1 node)
+  scc_groups <- scc_groups[sapply(scc_groups, length) > 1]
+
+  # 3. Assign colors
+  vertex_colors <- rep("grey80", vcount(g))
+  names(vertex_colors) <- V(g)$name
+
+  if (length(scc_groups) > 0) {
+
+    palette_colors <- rainbow(length(scc_groups))
+
+    for (i in seq_along(scc_groups)) {
+      vertex_colors[scc_groups[[i]]] <- palette_colors[i]
+    }
+  }
+
+  # 4. Plot
+  plot(g,
+       layout = layout_fun(g),
+       vertex.color = vertex_colors,
+       vertex.size = vertex_size,
+       vertex.label.cex = 0.8,
+       vertex.frame.color = "black",
+       edge.arrow.size = edge_arrow_size,
+       main = paste0("Enhancement recruitment loops (Nurse -> Recruit)"))
+
+  # Return invisibly for further manipulation if needed
+  invisible(list(graph = g,
+                 scc = scc_groups))
+}
diff --git a/tests/testthat/test-RN_dims.R b/tests/testthat/test-RN_dims.R
new file mode 100644
index 0000000..b9ba033
--- /dev/null
+++ b/tests/testthat/test-RN_dims.R
@@ -0,0 +1,98 @@
+test_that("RN_dims works for recruitment networks in Amoladeras", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  out <- RN_dims(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec"
+  )
+
+  # estructura del output
+  expect_s3_class(out, "data.frame")
+  expect_equal(rownames(out),
+               c("Num. Nodes", "Num. Links", "Connectance"))
+  expect_true(all(out$Value >= 0))
+})
+
+#-------------------
+test_that("RN_dims computes connectance correctly for rec networks", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  out <- RN_dims(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec"
+  )
+
+  # reconstrucción manual del cálculo
+  df <- comm_to_RN_UNI(Amoladeras_com, Amoladeras_cov)
+  n_nodes <- length(unique(c(df$Canopy, df$Recruit)))
+  n_links <- sum(df$Pcr)
+  connectance_expected <- n_links / (n_nodes^2 - n_nodes)
+
+  expect_equal(out["Connectance", "Value"], connectance_expected)
+})
+
+#---------------------
+test_that("RN_dims works for facilitation networks in Amoladeras", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  out <- RN_dims(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "fac"
+  )
+
+  expect_s3_class(out, "data.frame")
+  expect_equal(rownames(out),
+               c("Num. Nurse sp","Num. Facilitated sp", "Num. Links", "Connectance"))
+  expect_true(all(out$Value >= 0))
+})
+
+#---------------------
+test_that("RN_dims works for competition networks in Amoladeras", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+ out <- RN_dims(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "comp"
+  )
+
+  expect_s3_class(out, "data.frame")
+  expect_equal(rownames(out),
+               c("Num. Canopy depressing sp","Num. Recruit depressed sp", "Num. Links", "Connectance"))
+  expect_true(all(out$Value >= 0))
+})
+
+#-----------------------------
+test_that("RN_dims connectance for bipartite networks is correctly calculated", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  fac_dims <- RN_dims(Amoladeras_com, Amoladeras_cov, int_type = "fac")
+
+  mat <- RN_to_matrix(Amoladeras_com, Amoladeras_cov,
+                      int_type = "fac", weight = "Pcr")
+
+  n_links <- sum(mat)
+  connectance_expected <- n_links / (nrow(mat) * ncol(mat))
+
+  expect_equal(fac_dims["Connectance", "Value"], connectance_expected)
+})
diff --git a/tests/testthat/test-RN_heatmap.R b/tests/testthat/test-RN_heatmap.R
new file mode 100644
index 0000000..2286978
--- /dev/null
+++ b/tests/testthat/test-RN_heatmap.R
@@ -0,0 +1,14 @@
+#################################################
+#Tests for function: RN_heatmap
+#################################################
+
+test_that("RN_heatmat returns a ggplot object", {
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+  heatm_rec_Fcr <- RN_heatmap(Amoladeras_com,Amoladeras_cov, int_type="rec", weight ="Fcr")
+  expect_is(heatm_rec_Fcr, "ggplot")
+  heatm_fac_Ns <- RN_heatmap(Amoladeras_com,Amoladeras_cov, int_type="fac", weight ="Ns")
+  expect_is(heatm_fac_Ns, "ggplot")
+  heatm_comp_RII <- RN_heatmap(Amoladeras_com,Amoladeras_cov, int_type="comp", weight ="RII")
+  expect_is(heatm_comp_RII, "ggplot")
+})
\ No newline at end of file
diff --git a/tests/testthat/test-RN_to_matrix.R b/tests/testthat/test-RN_to_matrix.R
new file mode 100644
index 0000000..ce92623
--- /dev/null
+++ b/tests/testthat/test-RN_to_matrix.R
@@ -0,0 +1,188 @@
+library(testthat)
+
+
+test_that("RN_to_matrix works (Amoladeras)", {
+
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  net <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec",
+    weight     = "Fcr"
+  )
+
+  expected <- structure(
+    c(3L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 9L, 0L, 0L, 0L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 1L, 1L, 
+      0L, 4L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 1L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      1L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 
+      0L, 14L, 0L, 5L, 0L, 0L, 1L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 
+      3L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 15L, 0L, 4L, 0L, 7L, 0L, 0L, 
+      12L, 14L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
+      0L, 0L, 1L, 0L, 8L, 0L, 0L, 0L, 5L, 2L, 2L, 1L, 1L, 0L, 0L, 0L, 
+      5L, 7L, 5L, 35L, 4L, 9L, 0L, 6L, 1L, 1L, 113L, 0L, 168L, 2L, 
+      8L, 0L, 4L, 14L, 20L, 14L, 36L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 1L, 2L, 20L, 17L, 0L, 3L, 0L, 4L, 0L, 0L, 9L, 0L, 
+      42L, 0L, 12L, 0L, 0L, 7L, 4L, 0L, 0L, 0L, 23L, 3L, 4L, 7L, 0L, 
+      4L, 0L, 6L, 0L, 2L, 10L, 1L, 10L, 0L, 41L, 0L, 17L, 0L, 6L, 17L, 
+      13L, 12L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 
+      0L, 3L, 0L, 7L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 0L, 1L, 0L, 0L, 
+      0L, 0L, 2L, 2L, 19L, 0L, 2L, 0L, 0L, 36L, 0L, 3L, 0L, 3L, 0L, 
+      0L, 4L, 1L, 0L, 0L, 0L, 54L, 1L, 7L, 18L, 11L, 796L, 34L, 89L, 
+      6L, 20L, 80L, 0L, 1839L, 0L, 528L, 19L, 60L, 1L, 16L, 308L, 170L, 
+      1030L, 5L, 4L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 3L, 0L, 1L, 0L, 0L, 
+      5L, 0L, 3L, 0L, 5L, 0L, 0L, 8L, 9L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 5L, 0L, 1L, 0L, 0L, 0L, 0L, 4L, 0L, 0L, 0L, 0L, 4L, 
+      0L, 0L, 0L, 0L, 2L, 0L, 0L, 2L, 4L, 32L, 1L, 9L, 0L, 5L, 1L, 
+      0L, 6L, 0L, 46L, 0L, 0L, 0L, 0L, 17L, 5L, 14L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
+      0L, 0L, 0L, 6L, 3L, 0L, 0L, 0L, 0L, 0L, 8L, 0L, 7L, 2L, 3L, 0L, 
+      0L, 0L, 12L, 3L, 0L, 0L, 5L, 0L, 0L, 5L, 5L, 58L, 2L, 7L, 1L, 
+      1L, 2L, 0L, 43L, 0L, 138L, 0L, 11L, 0L, 5L, 54L, 19L, 34L, 1L, 
+      0L, 2L, 2L, 3L, 2L, 1L, 86L, 3L, 4L, 1L, 2L, 6L, 0L, 7L, 0L, 
+      133L, 4L, 27L, 0L, 0L, 87L, 65L, 9L, 0L, 0L, 0L, 0L, 0L, 4L, 
+      0L, 0L, 0L, 8L, 0L, 3L, 0L, 1L, 4L, 0L, 10L, 0L, 0L, 0L, 0L, 
+      3L, 2L, 5L, 1L, 0L, 0L, 0L, 4L, 0L, 32L, 12L, 4L, 4L, 0L, 13L, 
+      2L, 0L, 9L, 0L, 40L, 0L, 14L, 0L, 7L, 9L, 3L, 9L, 2L, 0L),
+    dim = c(24L, 24L),
+    dimnames = list(
+      c("Artemisia_barrelieri","Artemisia_campestris","Asparagus_albus",
+        "Asparagus_horridus","Ballota_hirsuta","Helichrysum_stoechas",
+        "Hyparrhenia_hirta","Launaea_arborescens","Launaea_lanifera",
+        "Lycium_intrincatum","Lygeum_spartum","Maytenus_senegalensis",
+        "Ononis_natrix","Open","Phagnalon_saxatile","Piptatherum_miliaceum",
+        "Salsola_oppositifolia","Stipa_tenacissima","Teucrium_lusitanicum",
+        "Teucrium_polium","Thymelaea_hirsuta","Thymus_hyemalis",
+        "Whitania_frutescens","Ziziphus_lotus"),
+      c("Artemisia_barrelieri","Artemisia_campestris","Asparagus_albus",
+        "Asparagus_horridus","Ballota_hirsuta","Helichrysum_stoechas",
+        "Hyparrhenia_hirta","Launaea_arborescens","Launaea_lanifera",
+        "Lycium_intrincatum","Lygeum_spartum","Maytenus_senegalensis",
+        "Ononis_natrix","Open","Phagnalon_saxatile","Piptatherum_miliaceum",
+        "Salsola_oppositifolia","Stipa_tenacissima","Teucrium_lusitanicum",
+        "Teucrium_polium","Thymelaea_hirsuta","Thymus_hyemalis",
+        "Whitania_frutescens","Ziziphus_lotus")
+    )
+  )
+
+  expect_equal(net, expected)
+})
+
+#-------------------------------
+test_that("RN_to_matrix rec has non-negative values for Fcr (Amoladeras)", {
+
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  net <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec",
+    weight     = "Fcr"
+  )
+
+  expect_true(all(net >= 0, na.rm = TRUE))
+})
+
+#-------------------------------
+
+test_that("RN_to_matrix fac produces sparse matrix (Amoladeras)", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  net <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "fac",
+    weight     = "Fcr"
+  )
+
+  # facilitation debería generar muchos ceros
+  prop_zero <- sum(net == 0) / length(net)
+
+  expect_true(prop_zero > 0.3)
+})
+
+#-------------------------------
+
+test_that("RN_to_matrix comp produces values only for depressing interactions", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  net <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "comp",
+    weight     = "Fcr"
+  )
+
+  # matriz con señal (no todo cero)
+  expect_true(sum(net != 0) > 0)
+})
+
+
+#-------------------------------
+
+test_that("Changing weight changes matrix values (Amoladeras)", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+
+  net_Fcr <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec",
+    weight     = "Fcr"
+  )
+
+  net_RII <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec",
+    weight     = "RII"
+  )
+
+  # misma dimensión
+  expect_equal(dim(net_Fcr), dim(net_RII))
+
+  # pero valores distintos
+  expect_false(identical(net_Fcr, net_RII))
+})
+
+#-------------------------------
+
+test_that("RN_to_matrix keeps Open column but it may be zero for some indices", {
+ 
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  net <- RN_to_matrix(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec",
+    weight     = "Ns"
+  )
+
+  expect_true("Open" %in% colnames(net))
+})
diff --git a/tests/testthat/test-associndex.R b/tests/testthat/test-associndex.R
new file mode 100644
index 0000000..847d5da
--- /dev/null
+++ b/tests/testthat/test-associndex.R
@@ -0,0 +1,82 @@
+# Test 1.- La función da error si no existe open
+test_that("associndex falla si no existe Open", {
+  int_no_open <- vent.int[vent.int$Canopy != "Open", ]
+  expect_error(
+    associndex(int_no_open, vent.cov),
+    "does not contain a node named Open"
+  )
+})
+# Test 2.- Devuelve null en la rama no allsp
+test_that("associndex devuelve NULL en combinación unused", {
+  expect_warning(
+    res <- associndex(vent.int, vent.cov,
+                      expand = "no",
+                      rm_sp_no_cover = "allsp")
+  )
+  expect_null(res)
+})
+
+#Test 3.- rama yes sp devuelve dataframe
+test_that("associndex modo rec devuelve dataframe", {
+  
+  res <- associndex(vent.int, vent.cov,
+                    expand = "yes",
+                    rm_sp_no_cover = "allsp")
+  
+  expect_s3_class(res, "data.frame")
+  expect_gt(nrow(res), 0)
+  
+})
+
+# Test 4.- Rama no onlycanopy devuelve dataframe
+test_that("associndex modo fac devuelve dataframe", {
+  
+  res <- associndex(vent.int, vent.cov,
+                    expand = "no",
+                    rm_sp_no_cover = "onlycanopy")
+  
+  expect_s3_class(res, "data.frame")
+  expect_gt(nrow(res), 0)
+  
+})
+
+# Test 5.- Rama yes onlycanopy devuelve dataframe
+test_that("associndex modo comp devuelve dataframe", {
+  
+  res <- associndex(vent.int, vent.cov,
+                    expand = "yes",
+                    rm_sp_no_cover = "onlycanopy")
+  
+  expect_s3_class(res, "data.frame")
+  expect_gt(nrow(res), 0)
+  
+})
+
+# Test 6.- la función calcula el thresold automáticamente
+test_that("associndex calcula threshold automáticamente", {
+  
+  expect_message(
+    associndex(vent.int, vent.cov,
+               expand = "yes",
+               rm_sp_no_cover = "allsp",
+               threshold_density = NULL),
+    "threshold density has been set"
+  )
+  
+})
+
+# Test 7.- La función no calcula el thresold si se da
+
+test_that("associndex no calcula threshold si se proporciona", {
+  
+  expect_no_message(
+    associndex(vent.int, vent.cov,
+               expand = "yes",
+               rm_sp_no_cover = "allsp",
+               threshold_density = 1000)
+  )
+  
+})
+
+
+
diff --git a/tests/testthat/test-canopy_service_test.R b/tests/testthat/test-canopy_service_test.R
new file mode 100644
index 0000000..ab257dd
--- /dev/null
+++ b/tests/testthat/test-canopy_service_test.R
@@ -0,0 +1,70 @@
+#################################################
+#Tests para canopy_service_test
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("canopy_service_test devuelve un data.frame con las columnas correctas", {
+  
+  res <- canopy_service_test(mysite_com, mysite_cov)
+  
+  expect_s3_class(res, "data.frame")
+  
+  expect_equal(
+    colnames(res),
+    c("Canopy", "Fc", "Ac", "Fro", "Ao",
+      "testability", "Significance", "Test_type", "Canopy_effect")
+  )
+  
+  expect_equal(nrow(res), 23)
+})
+
+
+#------------------------------------
+
+test_that("Valores agregados correctos para Artemisia_barrelieri", {
+  
+  res <- canopy_service_test(mysite_com, mysite_cov)
+  
+  fila <- res[res$Canopy == "Artemisia_barrelieri", ]
+  
+  expect_equal(fila$Fc, 18)
+  expect_equal(fila$Fc, sum(mysite_com[mysite_com$Canopy=="Artemisia_barrelieri" ,"Frequency"]))
+  expect_equal(fila$Fro, 5111)
+  expect_equal(fila$Fro, sum(mysite_com[mysite_com$Canopy=="Open" ,"Frequency"]))
+  expect_equal(fila$Ac, 1,tolerance = 1e-6)
+  expect_equal(fila$Ao, 6927,tolerance = 1e-6)
+
+
+})
+
+
+#--------------------------------------
+
+test_that("Clasificación Canopy_effect correcta para especies conocidas", {
+  
+  res <- canopy_service_test(mysite_com, mysite_cov)
+  
+  expect_equal(
+    res$Canopy_effect[res$Canopy == "Artemisia_barrelieri"],
+    "Facilitative"
+  )
+  
+  expect_equal(
+    res$Canopy_effect[res$Canopy == "Artemisia_campestris"],
+    "Neutral"
+  )
+})
+
diff --git a/tests/testthat/test-comm_to_RN.R b/tests/testthat/test-comm_to_RN.R
new file mode 100644
index 0000000..2c9d335
--- /dev/null
+++ b/tests/testthat/test-comm_to_RN.R
@@ -0,0 +1,21 @@
+#################################################
+#Tests for function: comm_to_RN
+#################################################
+
+test_that("comm_to_RN works", {
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+  comm_to_RN_rec <- comm_to_RN(Amoladeras_com, Amoladeras_cov, expand="yes",rm_sp_no_cover="allsp")
+  comm_to_RN_fac <- comm_to_RN(Amoladeras_com, Amoladeras_cov, expand="no",rm_sp_no_cover="onlycanopy")
+  comm_to_RN_comp <- comm_to_RN(Amoladeras_com, Amoladeras_cov, expand="yes",rm_sp_no_cover="onlycanopy")
+  
+  expect_equal(unlist(comm_to_RN_rec[1, 1:7]), c("Artemisia_barrelieri", "Artemisia_barrelieri", 3, 1, 1, 1, 1), check.attributes = FALSE)
+  expect_equal(dim(comm_to_RN_rec), c(576,7))
+  expect_equal(unlist(comm_to_RN_fac[1, 1:6]), c("Artemisia_barrelieri", "Artemisia_barrelieri", 3, 1, 1, 1), check.attributes = FALSE)
+  expect_equal(dim(comm_to_RN_fac), c(229,6))
+  expect_equal(unlist(comm_to_RN_comp[1, 1:6]), c("Artemisia_barrelieri", "Artemisia_barrelieri", 3, 1, 1, 1), check.attributes = FALSE)
+  expect_equal(dim(comm_to_RN_comp), c(624,6))
+  
+  expect_warning(comm_to_RN(Amoladeras_com, Amoladeras_cov, expand="no",rm_sp_no_cover="allsp"))
+  
+})
\ No newline at end of file
diff --git a/tests/testthat/test-cum_values.R b/tests/testthat/test-cum_values.R
new file mode 100644
index 0000000..55a8054
--- /dev/null
+++ b/tests/testthat/test-cum_values.R
@@ -0,0 +1,15 @@
+#################################################
+#Tests for function: cum_values
+#################################################
+
+test_that("cum_values works", {
+  Amoladeras <- comm_subset(RecruitNet, site = "Amoladeras")
+  cum_ecount <- cum_values(Amoladeras, property="ecount", k=10)
+  cum_vcount <- cum_values(Amoladeras, property="vcount", k=10)
+  cum_edge_density <- cum_values(Amoladeras, property="edge_density", k=10)
+  expect_equal(dim(cum_ecount$Data), c(200,2))
+  expect_equal(length(unique(cum_ecount$Data$sampleSize)),20)
+  expect_equal(cum_ecount$Data$Value[200], 229)
+  expect_equal(cum_vcount$Data$Value[200], 26)
+  expect_equal(cum_edge_density$Data$Value[200], 0.3523077)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-funtopol_rec.R b/tests/testthat/test-funtopol_rec.R
new file mode 100644
index 0000000..544a89b
--- /dev/null
+++ b/tests/testthat/test-funtopol_rec.R
@@ -0,0 +1,82 @@
+#################################################
+#Tests para funtopol_rec
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("funtopol_rec devuelve la estructura correcta", {
+  
+  res <- funtopol_rec(mysite_com, mysite_cov)
+  
+  expect_type(res, "list")
+  expect_named(res, c("Descriptors", "Functional_classification"))
+  
+  expect_s3_class(res$Descriptors, "data.frame")
+  expect_s3_class(res$Functional_classification, "data.frame")
+})
+
+#--------------------------------------------------
+
+test_that("funtopol_rec devuelve los valores esperados en Descriptors", {
+  
+  res <- funtopol_rec(mysite_com, mysite_cov)
+  df  <- res$Descriptors
+  
+  expect_equal(df["Num. nodes", "Value"], 24)
+  expect_equal(df["Num. edges", "Value"], 221)
+  expect_equal(df["Connectance", "Value"], 0.4, tolerance = 1e-6)
+  expect_equal(df["Num. non-trivial SCCs", "Value"], 1)
+  expect_equal(df["Num. core species", "Value"], 19)
+})
+
+#--------------------------------------------------
+
+test_that("Clasificación funcional tiene el número correcto de especies", {
+  
+  res <- funtopol_rec(mysite_com, mysite_cov)
+  fc  <- res$Functional_classification
+  
+  a<-sort(a$Functional_classification$id)
+  b<-sort(intersect(unique(mysite_com$Recruit),unique(mysite_cov$Canopy)))
+  expect_equal(nrow(fc), 23)
+  expect_true(all(c("id", "group") %in% colnames(fc)))
+  expect_equal(a, b)
+
+})
+
+##--------------------------------------------------
+
+test_that("Artemisia_barrelieri está clasificada como Core", {
+  
+  res <- funtopol_rec(mysite_com, mysite_cov)
+  fc  <- res$Functional_classification
+  
+  fila <- fc[fc$id == "Artemisia_campestris", ]
+  
+  expect_equal(fila$group, "Satellite")
+})
+
+# --------------------------------------------------
+
+test_that("Da error si no existe Open", {
+  
+  com_sin_open <- mysite_com
+  com_sin_open$Canopy <- gsub("Open", "OPEN_WRONG", com_sin_open$Canopy)
+  
+  expect_error(
+    funtopol_rec(com_sin_open, mysite_cov),
+    "does not contain a node named Open"
+  )
+})
diff --git a/tests/testthat/test-int_significance.R b/tests/testthat/test-int_significance.R
new file mode 100644
index 0000000..0fad6b8
--- /dev/null
+++ b/tests/testthat/test-int_significance.R
@@ -0,0 +1,98 @@
+library(testthat)
+
+data(RecruitNet)
+data(CanopyCover)
+
+test_that("int_significance works for recruitment network in Amoladeras", {
+  
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  out <- int_significance(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "rec"
+  )
+
+  expect_s3_class(out, "data.frame")
+  expect_true(nrow(out) > 0)
+  expect_true("Test_type" %in% colnames(out))
+})
+
+#----------------------
+test_that("int_significance works for facilitation network in Amoladeras", {
+ #data(RecruitNet)
+ #data(CanopyCover)
+
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  out <- int_significance(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "fac"
+  )
+
+  expect_s3_class(out, "data.frame")
+  expect_true(all(out$Effect_int == "Enhancing"))
+})
+
+#----------------
+test_that("int_significance works for competition network in Amoladeras", {
+  
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  out <- int_significance(
+    int_data   = Amoladeras_com,
+    cover_data = Amoladeras_cov,
+    int_type   = "comp"
+  )
+
+  expect_s3_class(out, "data.frame")
+  expect_true(all(out$Effect_int == "Depressing"))
+})
+
+#----------------
+
+test_that("int_significance fails in Amoladeras if 'Open' canopy is missing", {
+  
+  #data(RecruitNet)
+  #data(CanopyCover)
+
+   Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+   Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+  int_data <- Amoladeras_com
+  int_data$Canopy <- as.character(int_data$Canopy)
+  int_data <- int_data[int_data$Canopy != "Open", ]
+
+  expect_error(
+    int_significance(
+      int_data   = int_data,
+      cover_data = Amoladeras_cov,
+      int_type   = "rec"
+    ),
+    "ERROR: tests cannot be conducted"
+  )
+})
+
+#--------------------
+test_that("int_significance message appears when multiple statistical tests are used (Amoladeras)", {
+  #data(RecruitNet)
+  #data(CanopyCover)
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+
+
+  expect_message(
+    int_significance(
+      int_data   =   Amoladeras_com,
+      cover_data =   Amoladeras_cov,
+      int_type   = "rec"
+    ),
+    "Different tests were used",
+    all = FALSE
+  )
+})
diff --git a/tests/testthat/test-node_degrees.R b/tests/testthat/test-node_degrees.R
new file mode 100644
index 0000000..96c121f
--- /dev/null
+++ b/tests/testthat/test-node_degrees.R
@@ -0,0 +1,23 @@
+#################################################
+#Tests for function: node_degrees
+#################################################
+
+test_that("node_degrees works", {
+  Amoladeras_com <- comm_subset(RecruitNet, site = "Amoladeras")
+  Amoladeras_cov <- comm_subset(CanopyCover, site = "Amoladeras")
+  node_degree_rec <- node_degrees(Amoladeras_com, Amoladeras_cov, int_type = "rec")
+  node_degree_fac <- node_degrees(Amoladeras_com, Amoladeras_cov, int_type = "fac")
+  node_degree_comp <- node_degrees(Amoladeras_com, Amoladeras_cov, int_type = "comp")
+  
+  expect_equal(unlist(node_degree_rec[1, 1:8]), c("Artemisia_barrelieri", 1, 5, 18, 3.85179402240326, 7, 90, 3.15278474320294), check.attributes = FALSE)
+  expect_equal(dim(node_degree_rec), c(24,8))
+  expect_equal(unlist(node_degree_fac$Canopy[1, 1:3]), c("Artemisia_barrelieri", 1, 4), check.attributes = FALSE)
+  expect_equal(dim(node_degree_fac$Canopy), c(18,3))
+  expect_equal(unlist(node_degree_fac$Recruit[1, 1:3]), c("Artemisia_barrelieri", 90, 3), check.attributes = FALSE)
+  expect_equal(dim(node_degree_fac$Recruit), c(19,3))
+  expect_equal(unlist(node_degree_comp$Canopy[1, 1:3]), c("Asparagus_albus", 27.725, 2), check.attributes = FALSE)
+  expect_equal(dim(node_degree_comp$Canopy), c(13,3))
+  expect_equal(unlist(node_degree_comp$Recruit[1, 1:3]), c("Artemisia_barrelieri", 90, 2), check.attributes = FALSE)
+  expect_equal(dim(node_degree_comp$Recruit), c(7,3))
+  
+})
\ No newline at end of file
diff --git a/tests/testthat/test-node_topol.R b/tests/testthat/test-node_topol.R
new file mode 100644
index 0000000..9d5a14f
--- /dev/null
+++ b/tests/testthat/test-node_topol.R
@@ -0,0 +1,126 @@
+#################################################
+#Tests para node_topol
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("node_topol devuelve NULL si int_type no es válido", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "wrong")
+  expect_null(res)
+})
+
+#------------------------------------
+
+test_that("node_topol rec devuelve estructura correcta", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "rec")
+  
+  expect_s3_class(res, "data.frame")
+  expect_equal(nrow(res), 24)
+  
+  expect_equal(
+    colnames(res),
+    c("Eigenvector centrality",
+      "Extended canopy service",
+      "Extended recruitment niche")
+  )
+})
+
+
+#------------------------------------
+test_that("node_topol fac devuelve estructura correcta", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "fac")
+  
+  expect_s3_class(res, "data.frame")
+  expect_equal(nrow(res), 21)
+  
+  expect_equal(
+    colnames(res),
+    c("Eigenvector nurse centrality",
+      "Extended nurse service",
+      "Extended facilitated niche")
+  )
+})
+
+#------------------------------------
+
+test_that("node_topol comp devuelve estructura correcta", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "comp")
+  
+  expect_s3_class(res, "data.frame")
+  expect_equal(nrow(res), 17)
+  
+  expect_equal(
+    colnames(res),
+    c("Eigenvector canopy centrality",
+      "Extended canopy depression effect",
+      "Extended recruitment depression ")
+  )
+})
+
+
+#------------------------------------
+
+test_that("Valores agregados correctos para Artemisia_barrelieri en rec", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "rec")
+
+  
+  fila <- res[rownames(res) == "Artemisia_barrelieri", ]
+  
+  expect_equal(fila$Eigenvector centrality, 0.384,tolerance = 1e-6)
+  expect_equal(fila$Extended canopy service, 20)
+  expect_equal(fila$Extended recruitment niche, 20)
+
+})
+
+#------------------------------------
+
+test_that("Valores agregados correctos para Artemisia_barrelieri en fac", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "fac")
+
+  
+  fila <- res[rownames(res) == "Artemisia_barrelieri", ]
+  
+  expect_equal(fila$Eigenvector centrality, 0.3921,tolerance = 1e-6)
+  expect_equal(fila$Extended canopy service, 18)
+  expect_equal(fila$Extended recruitment niche, 17)
+
+})
+
+#------------------------------------
+
+test_that("Valores agregados correctos para Artemisia_barrelieri en comp", {
+  
+  res <- node_topol(mysite_com, mysite_cov, int_type = "comp")
+
+  
+  fila <- res[rownames(res) == "Artemisia_barrelieri", ]
+  
+  expect_equal(fila$Eigenvector centrality, 0,tolerance = 1e-6)
+  expect_equal(fila$Extended canopy service, 0)
+  expect_equal(fila$Extended recruitment niche, 2)
+
+})
+
+
+
+
+
+
diff --git a/tests/testthat/test-recruitment_niche_test.R b/tests/testthat/test-recruitment_niche_test.R
new file mode 100644
index 0000000..b924c60
--- /dev/null
+++ b/tests/testthat/test-recruitment_niche_test.R
@@ -0,0 +1,61 @@
+# Test 1.- La función devuelve un dataframe válido
+test_that("recruitment_niche_test devuelve dataframe válido", {
+  res <- recruitment_niche_test(vent.int, vent.cov)
+  expect_s3_class(res, "data.frame")
+  expect_true(all(c(
+    "Recruit",
+    "Fr",
+    "Av",
+    "Fro",
+    "Ao",
+    "testability",
+    "Significance",
+    "Test_type",
+    "Veg_effect"
+  ) %in% names(res)))
+})
+
+# Test 2.- Test type y Veg effect contienen categorías válidas
+test_that("Test_type y Veg_effect contienen categorías válidas", {
+  res <- recruitment_niche_test(vent.int, vent.cov)
+  expect_true(all(res$Test_type %in% c("Binomial", "Chi-square")))
+  expect_true(all(res$Veg_effect %in% c(
+    
+    "Facilitated",
+    "Depressed",
+    "Neutral",
+    "Not testable"
+  )))
+  expect_false(any(is.na(res$Test_type)))
+  expect_false(any(is.na(res$Veg_effect)))
+})
+
+# Test 3.- coherencia testability y veg_effect
+test_that("Coherencia completa entre testability y Veg_effect", {
+  res <- recruitment_niche_test(vent.int, vent.cov)
+  # Not testable
+  idx_nt <- which(res$testability > 0.05)
+  if(length(idx_nt) > 0){
+    expect_true(all(res$Veg_effect[idx_nt] == "Not testable"))
+  }
+  # Facilitated
+  idx_fac <- which(res$Veg_effect == "Facilitated")
+  if(length(idx_fac) > 0){
+    expect_true(all((res$Fr[idx_fac]/res$Av[idx_fac]) >
+                      (res$Fro[idx_fac]/res$Ao[idx_fac])))
+  }
+  # Depressed
+  idx_dep <- which(res$Veg_effect == "Depressed")
+  if(length(idx_dep) > 0){
+    expect_true(all((res$Fr[idx_dep]/res$Av[idx_dep]) <
+                      (res$Fro[idx_dep]/res$Ao[idx_dep])))
+  }
+})
+
+# Test 4.- se lanza el warning sobre diferentes tipos de tests  
+test_that("lanza warning cuando hay varios tipos de tests", {
+  expect_warning(
+    recruitment_niche_test(vent.int, vent.cov),
+    "Different tests were used")
+})
+
diff --git a/tests/testthat/test-topol_depre.R b/tests/testthat/test-topol_depre.R
new file mode 100644
index 0000000..e79b759
--- /dev/null
+++ b/tests/testthat/test-topol_depre.R
@@ -0,0 +1,126 @@
+#################################################
+#Tests para topol_depre
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("topol_depre returns correct list structure (real data)", {
+  
+  res <- topol_depre(mysite_com, mysite_cov, direction = "out")
+  
+  expect_type(res, "list")
+  expect_named(res, c("loops", "simple"))
+  expect_named(res$loops, c("summary", "nodes"))
+  expect_named(res$simple, c("summary", "nodes"))
+})
+
+#------------------------------------
+
+test_that("Real dataset (Amoladeras) contains no reciprocal loops", {
+  
+  res <- topol_depre(mysite_com, mysite_cov, direction = "out")
+  
+  expect_equal(nrow(res$loops$summary), 0)
+  expect_length(res$loops$nodes, 0)
+})
+
+#------------------------------------
+
+test_that("Dummy dataset detects strongly connected components (loops)", {
+  
+  res <- topol_depre(test_data$com, test_data$cov, direction = "out")
+  
+  expect_true(nrow(res$loops$summary) > 0)
+  expect_true(length(res$loops$nodes) > 0)
+  
+  # Each SCC should contain more than one node
+  expect_true(all(res$loops$summary$n_nodos > 1))
+})
+
+#------------------------------------
+
+test_that("Simple paths (out) have correct structure", {
+  
+  res <- topol_depre(mysite_com, mysite_cov, direction = "out")
+  
+  expect_true(nrow(res$simple$summary) > 0)
+  expect_true(all(res$simple$summary$n_nodes_in_path >= 1))
+  
+  for (node in names(res$simple$nodes)) {
+    paths <- res$simple$nodes[[node]]
+    for (p in paths) {
+      expect_equal(p[1], node)
+    }
+  }
+})
+
+#------------------------------------
+
+test_that("Simple paths (in) have correct structure", {
+  
+  res <- topol_depre(mysite_com, mysite_cov, direction = "in")
+  
+  expect_type(res, "list")
+  expect_true(nrow(res$simple$summary) > 0)
+})
+
+#------------------------------------
+
+test_that("Dummy dataset works for both directions", {
+  
+  res_out <- topol_depre(test_data$com, test_data$cov, direction = "out")
+  res_in  <- topol_depre(test_data$com, test_data$cov, direction = "in")
+  
+  expect_true(nrow(res_out$simple$summary) > 0)
+  expect_true(nrow(res_in$simple$summary) > 0)
+})
+
+#------------------------------------
+
+test_that("No path is a subpath of another (out direction)", {
+  
+  res <- topol_depre(mysite_com, mysite_cov, direction = "out")
+  
+  es_subpath <- function(p, q) {
+    lp <- length(p)
+    lq <- length(q)
+    if (lp >= lq) return(FALSE)
+    for (i in 1:(lq - lp + 1)) {
+      if (all(q[i:(i + lp - 1)] == p)) return(TRUE)
+    }
+    FALSE
+  }
+  
+  for (paths in res$simple$nodes) {
+    if (length(paths) > 1) {
+      for (i in seq_along(paths)) {
+        for (j in seq_along(paths)) {
+          if (i != j) {
+            expect_false(es_subpath(paths[[i]], paths[[j]]))
+          }
+        }
+      }
+    }
+  }
+})
+
+#------------------------------------
+
+test_that("Invalid direction argument throws error", {
+  
+  expect_error(
+    topol_depre(mysite_com, mysite_cov, direction = "sideways")
+  )
+})
diff --git a/tests/testthat/test-topol_fac.R b/tests/testthat/test-topol_fac.R
new file mode 100644
index 0000000..9f922ba
--- /dev/null
+++ b/tests/testthat/test-topol_fac.R
@@ -0,0 +1,119 @@
+#################################################
+#Tests para topol_fac
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("topol_fac provides teh correct structure", {
+  res <- topol_fac(mysite_com, mysite_cov, direction = "out")
+  
+  expect_type(res, "list")
+  expect_named(res, c("loops", "simple"))
+  
+  # loops
+  expect_named(res$loops, c("summary", "nodes"))
+  expect_s3_class(res$loops$summary, "data.frame")
+  expect_type(res$loops$nodes, "list")
+  
+  # simple
+  expect_named(res$simple, c("summary", "nodes"))
+  expect_s3_class(res$simple$summary, "data.frame")
+  expect_type(res$simple$nodes, "list")
+})
+
+#------------------------------------
+test_that("loops summary coincides with the actual nodes ", {
+  res <- topol_fac(mysite_com, mysite_cov, direction = "out")
+  
+  expect_equal(nrow(res$loops$summary), 1)
+  
+  scc_id <- res$loops$summary$scc_id
+  n_nodos_reportado <- res$loops$summary$n_nodos
+  
+  nodos_reales <- res$loops$nodes[[as.character(scc_id)]]
+  
+  expect_equal(length(nodos_reales), n_nodos_reportado)
+})
+
+#------------------------------------
+
+test_that("simple summary coincides with real paths", {
+  res <- topol_fac(mysite_com, mysite_cov, direction = "out")
+  
+  summary_df <- res$simple$summary
+  nodes_list <- res$simple$nodes
+  
+  for(i in seq_len(nrow(summary_df))) {
+    
+    nodo <- summary_df$nodo[i]
+    path_index <- summary_df$path_index[i]
+    n_nodes_in_path <- summary_df$n_nodes_in_path[i]
+    
+    path_real <- nodes_list[[nodo]][[path_index]]
+    
+    expect_equal(length(path_real), n_nodes_in_path)
+  }
+})
+
+#------------------------------------
+
+test_that("simple paths start by its original node", {
+  res <- topol_fac(mysite_com, mysite_cov, direction = "out")
+  
+  for(nodo in names(res$simple$nodes)) {
+    for(path in res$simple$nodes[[nodo]]) {
+      expect_equal(path[1], nodo)
+    }
+  }
+})
+
+#------------------------------------
+
+test_that("simple paths do not have repeted nodes", {
+  res <- topol_fac(mysite_com, mysite_cov, direction = "out")
+  
+  for(nodo in names(res$simple$nodes)) {
+    for(path in res$simple$nodes[[nodo]]) {
+      expect_equal(length(path), length(unique(path)))
+    }
+  }
+})
+
+#------------------------------------
+
+test_that("terminal species in a path has a length path of 1", {
+  res <- topol_fac(mysite_com, mysite_cov, direction="out")
+  
+  terminales <- c("Artemisia_campestris",
+                  "Maytenus_senegalensis",
+                  "Teucrium_lusitanicum")
+  
+  for(sp in terminales) {
+    expect_equal(
+      res$simple$summary$n_nodes_in_path[
+        res$simple$summary$nodo == sp
+      ],
+      1
+    )
+  }
+})
+
+#------------------------------------
+
+test_that("all nodes in summary are nodes", {
+  res <- topol_fac(mysite_com, mysite_cov, direction="out")
+  
+  expect_true(all(res$simple$summary$nodo %in% names(res$simple$nodes)))
+})
diff --git a/tests/testthat/test-visu_funtopol_rec.R b/tests/testthat/test-visu_funtopol_rec.R
new file mode 100644
index 0000000..f104752
--- /dev/null
+++ b/tests/testthat/test-visu_funtopol_rec.R
@@ -0,0 +1,97 @@
+#################################################
+#Tests para visu_funtopol_rec
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("Los nodos Core en la visualización coinciden con funtopol_rec()", {
+  
+  # Output ecológico
+  ft <- funtopol_rec(mysite_com, mysite_cov)$Functional_classification
+  core_expected <- ft$id[ft$group == "Core"]
+  
+  # Visualización
+  net <- visu_funtopol_rec(mysite_com, mysite_cov)
+  nodes_vis <- net$x$nodes
+  
+  core_visual <- nodes_vis$id[nodes_vis$group == "Core"]
+  
+  expect_setequal(core_visual, core_expected)
+})
+
+#------------------------------------
+
+test_that("Los nodos Core en la visualización coinciden con funtopol_rec()", {
+  
+  # Output ecológico
+  ft <- funtopol_rec(mysite_com, mysite_cov)$Functional_classification
+  core_expected <- ft$id[ft$group == "Satellite"]
+  
+  # Visualización
+  net <- visu_funtopol_rec(mysite_com, mysite_cov)
+  nodes_vis <- net$x$nodes
+  
+  core_visual <- nodes_vis$id[nodes_vis$group == "Satellite"]
+  
+  expect_setequal(core_visual, core_expected)
+})
+
+#------------------------------------
+
+test_that("El grupo Core está asignado al color verde correcto", {
+  
+  net <- visu_funtopol_rec(mysite_com, mysite_cov)
+  
+  groups <- net$x$options$groups
+  
+  expect_true("Core" %in% names(groups))
+  expect_equal(groups$Open$color$background, "#F0E442")
+  expect_equal(groups$Core$color$background, "#009E73")
+  expect_equal(groups$Satellite$color$background, "#0072B2")
+  expect_equal(groups$Strict_transients$color$background, "#D55E00")
+  expect_equal(groups$Disturbance_dependent_transients$color$background, "#CC79A7")
+})
+
+#------------------------------------
+
+test_that("No hay nodos mal clasificados como Core en la visualización", {
+  
+  ft <- funtopol_rec(mysite_com, mysite_cov)$Functional_classification
+  non_core <- ft$id[ft$group != "Core"]
+  
+  net <- visu_funtopol_rec(mysite_com, mysite_cov)
+  nodes_vis <- net$x$nodes
+  
+  wrongly_core <- intersect(
+    nodes_vis$id[nodes_vis$group == "Core"],
+    non_core
+  )
+  
+  expect_length(wrongly_core, 0)
+})
+
+#------------------------------------
+
+test_that("Todos los nodos Core aparecen en la visualización", {
+  
+  core_expected <- funtopol_rec(mysite_com, mysite_cov)$Functional_classification |>
+    subset(group == "Core") |>
+    (\(x) x$id)()
+  
+  net <- visu_funtopol_rec(mysite_com, mysite_cov)
+  
+  expect_true(all(core_expected %in% net$x$nodes$id))
+})
+
diff --git a/tests/testthat/test-visu_net.R b/tests/testthat/test-visu_net.R
new file mode 100644
index 0000000..542b08a
--- /dev/null
+++ b/tests/testthat/test-visu_net.R
@@ -0,0 +1,51 @@
+
+# Test 1.- la función valida o no argumentos erróneos
+test_that("visu_net valida argumentos correctamente", {
+  expect_error(visu_net(vent.int, vent.cov,
+                        int_type = "cherrypickers"))
+  expect_error(visu_net(vent.int, vent.cov,
+                        weight = "machupichu"))
+  expect_error(visu_net(vent.int, vent.cov,
+                        mode = "cucurucho"))
+})
+
+# Test 2.- La función lanza el warning rec bi
+test_that("visu_net lanza warning para rec + bi", {
+  expect_warning(
+    visu_net(vent.int, vent.cov,
+             int_type = "rec",
+             mode = "bi"),
+    "recruitment networks"
+  )
+})
+
+# Test 3.- la función de detiene cunado tiene valores weight incompatibles
+test_that("visu_net detiene si weight incompatible en rec", {
+  expect_error(
+    visu_net(vent.int, vent.cov,
+             int_type = "rec",
+             weight = "Ns",
+             mode = "uni"),
+    "inconsistency"
+  )
+})
+
+# Test 4.-  La función  devuelve un objeto igraph unipartito
+test_that("visu_net devuelve objeto igraph en modo unipartito", {
+  skip_if_not_installed("igraph")
+  g <- visu_net(vent.int, vent.cov,
+                int_type = "fac",
+                mode = "uni")
+  expect_s3_class(g, "igraph")
+  # Opcional pero recomendable:
+  expect_true(igraph::vcount(g) > 0)
+  expect_true(igraph::ecount(g) > 0)
+})
+# Test 5.- La función no devuelve igraph, no tiene return() ==> espera un NULL
+test_that("visu_net modo bi no devuelve igraph", {
+  skip_if_not_installed("bipartite")
+  res <- visu_net(vent.int, vent.cov,
+                  int_type = "fac",
+                  mode = "bi")
+  expect_null(res)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-visu_topol_depre.R b/tests/testthat/test-visu_topol_depre.R
new file mode 100644
index 0000000..c19ae2c
--- /dev/null
+++ b/tests/testthat/test-visu_topol_depre.R
@@ -0,0 +1,105 @@
+#################################################
+#Tests para visu_topol_depre
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("visu_topol_depre returns graph and SCC list", {
+  
+  expect_silent({
+    res <- visu_topol_depre(test_data$com,
+                            test_data$cov,
+                            direction = "out")
+  })
+  
+  expect_type(res, "list")
+  expect_named(res, c("graph", "scc"))
+  expect_s3_class(res$graph, "igraph")
+})
+
+#------------------------------------
+
+test_that("Graph contains correct number of nodes", {
+  
+  res <- suppressWarnings(
+    visu_topol_depre(test_data$com,
+                     test_data$cov,
+                     direction = "out")
+  )
+  
+  expected_species <- unique(c(test_data$com$Canopy,
+                               test_data$com$Recruit))
+  
+  expect_equal(vcount(res$graph), length(expected_species))
+})
+
+#------------------------------------
+
+test_that("SCC detection works when loops exist", {
+  
+  res <- suppressWarnings(
+    visu_topol_depre(test_data$com,
+                     test_data$cov,
+                     direction = "out")
+  )
+  
+  # There should be at least one SCC with >1 node
+  expect_true(length(res$scc) > 0)
+  expect_true(all(sapply(res$scc, length) > 1))
+})
+
+#------------------------------------
+
+test_that("Works with real dataset even if no loops exist", {
+  
+  expect_silent({
+    res <- visu_topol_depre(mysite_com,
+                            mysite_cov,
+                            direction = "out")
+  })
+  
+  expect_type(res$scc, "list")
+})
+
+
+#------------------------------------
+
+
+test_that("Custom layout function does not break plotting", {
+  
+  expect_silent(
+    visu_topol_depre(test_data$com,
+                     test_data$cov,
+                     direction = "out",
+                     layout_fun = layout_in_circle)
+  )
+})
+
+
+#------------------------------------
+
+test_that("Custom vertex and edge parameters work", {
+  
+  expect_silent(
+    visu_topol_depre(test_data$com,
+                     test_data$cov,
+                     direction = "out",
+                     vertex_size = 10,
+                     edge_arrow_size = 0.2)
+  )
+})
+
+
+
diff --git a/tests/testthat/test-visu_topol_fac.R b/tests/testthat/test-visu_topol_fac.R
new file mode 100644
index 0000000..fcd402f
--- /dev/null
+++ b/tests/testthat/test-visu_topol_fac.R
@@ -0,0 +1,125 @@
+#################################################
+#Tests para visu_topol_fac
+#################################################
+
+#library(testthat)
+
+#load data
+mypath<-getwd()
+download_RN() # Run only the first time you use the package.
+setwd(mypath)
+RecruitNet <-read.csv("RecruitNet.csv")
+CanopyCover <-read.csv("CanopyCover.csv")
+
+mysite_com <- comm_subset_UNI(RecruitNet, "Amoladeras")
+mysite_cov <- comm_subset_UNI(CanopyCover, "Amoladeras")
+
+#------------------------------------
+
+test_that("visu_topol_fac returns graph and SCC list", {
+  
+  expect_silent({
+    res <- visu_topol_fac(test_data$com,
+                          test_data$cov,
+                          direction = "out")
+  })
+  
+  expect_type(res, "list")
+  expect_named(res, c("graph", "scc"))
+  expect_s3_class(res$graph, "igraph")
+  expect_type(res$scc, "list")
+})
+
+#------------------------------------
+
+test_that("Graph has correct number of vertices", {
+  
+  res <- suppressWarnings(
+    visu_topol_fac(test_data$com,
+                   test_data$cov,
+                   direction = "out")
+  )
+  
+  expected_species <- unique(c(test_data$com$Canopy,
+                               test_data$com$Recruit))
+  
+  expect_equal(vcount(res$graph), length(expected_species))
+})
+
+#------------------------------------
+
+test_that("SCC groups match those detected by fac_topol", {
+  
+  res_visu <- suppressWarnings(
+    visu_topol_fac(test_data$com,
+                   test_data$cov,
+                   direction = "out")
+  )
+  
+  res_topol <- fac_topol(test_data$com,
+                         test_data$cov,
+                         direction = "out")
+  
+  # Compare number of SCCs
+  expect_equal(length(res_visu$scc),
+               length(res_topol$loops$nodes))
+  
+  # Compare sizes of SCCs
+  visu_sizes <- sort(sapply(res_visu$scc, length))
+  topol_sizes <- sort(sapply(res_topol$loops$nodes, length))
+  
+  expect_equal(visu_sizes, topol_sizes)
+})
+#------------------------------------
+
+test_that("Detects SCCs with more than one node", {
+  
+  res <- suppressWarnings(
+    visu_topol_fac(test_data$com,
+                   test_data$cov,
+                   direction = "out")
+  )
+  
+  expect_true(length(res$scc) > 0)
+  expect_true(all(sapply(res$scc, length) > 1))
+})
+
+#------------------------------------
+
+test_that("Custom layout function works correctly", {
+  
+  expect_silent(
+    visu_topol_fac(test_data$com,
+                   test_data$cov,
+                   direction = "out",
+                   layout_fun = layout_in_circle)
+  )
+})
+
+#------------------------------------
+
+test_that("Custom vertex and edge parameters do not break plotting", {
+  
+  expect_silent(
+    visu_topol_fac(test_data$com,
+                   test_data$cov,
+                   direction = "out",
+                   vertex_size = 10,
+                   edge_arrow_size = 0.2)
+  )
+})
+
+
+
+#------------------------------------
+
+test_that("Works with real dataset without errors", {
+  
+  expect_silent({
+    res <- visu_topol_fac(mysite_com,
+                          mysite_cov,
+                          direction = "out")
+  })
+  
+  expect_s3_class(res$graph, "igraph")
+})