diff --git a/R/adjmat.r b/R/adjmat.r
index 06c4d31..c0b08c9 100644
--- a/R/adjmat.r
+++ b/R/adjmat.r
@@ -464,29 +464,59 @@ adjmat_to_edgelist.list <- function(graph, undirected, keep.isolates) {
 #' @author George G. Vega Yon & Thomas W. Valente
 toa_mat <- function(obj, labels=NULL, t0=NULL, t1=NULL) {
 
+  if (inherits(obj, "matrix")) {
+    num_of_behaviors <- dim(obj)[2]
+  } else if (inherits(obj, "diffnet")){
+    if (inherits(obj$toa, "matrix")) {
+      num_of_behaviors <- dim(obj$toa)[2]}
+    else {num_of_behaviors <- 1}
+  } else {num_of_behaviors <- 1}
+
   if (!inherits(obj, "diffnet")) {
     if (!length(t0)) t0 <- min(obj, na.rm = TRUE)
     if (!length(t1)) t1 <- max(obj, na.rm = TRUE)
   }
 
-  cls <- class(obj)
-  ans <- if ("numeric" %in% cls) {
-    toa_mat.numeric(obj, labels, t0, t1)
-    } else if ("integer" %in% cls) {
-    toa_mat.integer(obj, labels, t0, t1)
-    } else if  ("diffnet" %in% cls) {
-    with(obj, list(adopt=adopt,cumadopt=cumadopt))
-    } else
-      stopifnot_graph(obj)
-
-
-  if (inherits(obj, "diffnet")) {
-    dimnames(ans$adopt) <- with(obj$meta, list(ids,pers))
-    dimnames(ans$cumadopt) <- with(obj$meta, list(ids,pers))
+  ans <- list()
+  if (num_of_behaviors == 1) {
+    cls <- class(obj)
+    ans[[1]] <- if ("numeric" %in% cls) {
+            toa_mat.numeric(obj, labels, t0, t1)
+            } else if ("integer" %in% cls) {
+            toa_mat.integer(obj, labels, t0, t1)
+            } else if  ("diffnet" %in% cls) {
+            with(obj, list(adopt=adopt,cumadopt=cumadopt))
+            } else {
+              stopifnot_graph(obj)
+            }
+  } else {
+    for (q in 1:num_of_behaviors) {
+      ans[[q]] <- if ("matrix" %in% class(obj)) {
+              if ("integer" %in% class(obj[,q])){
+                toa_mat.integer(obj[,q], labels, t0, t1)
+              } else if ("numeric" %in% class(obj[,q])) { # Why included?
+                toa_mat.numeric(obj[,q], labels, t0, t1)
+              }
+            } else if  ("diffnet" %in% class(obj)) { # Why included?
+              with(obj, list(adopt=adopt[[q]],cumadopt=cumadopt[[q]]))
+            } else {
+              stopifnot_graph(obj[,q])
+            }
+    }
   }
 
+  for (q in 1:num_of_behaviors) {
+    if (inherits(obj, "diffnet")) {
+      dimnames(ans[[q]]$adopt) <- with(obj$meta, list(ids,pers))
+      dimnames(ans[[q]]$cumadopt) <- with(obj$meta, list(ids,pers))
+    }
+  }
 
-  return(ans)
+  if (num_of_behaviors==1) {
+    return(ans[[1]])
+  } else {
+    return(ans)
+  }
 }
 
 toa_mat.default <- function(per, t0, t1) {
diff --git a/R/diffnet-class.r b/R/diffnet-class.r
index 6704767..87163b3 100644
--- a/R/diffnet-class.r
+++ b/R/diffnet-class.r
@@ -556,6 +556,12 @@ new_diffnet <- function(graph, toa, t0=min(toa, na.rm = TRUE), t1=max(toa, na.rm
     return(graph)
   }
 
+  # Step 0.1: Setting num_of_behavior ------------------------------------------
+
+  if (inherits(toa, "matrix")) {
+    num_of_behaviors <- dim(toa)[2]
+  } else {num_of_behaviors <- 1}
+
   # Step 1.1: Check graph ------------------------------------------------------
   meta <- classify_graph(graph)
   if (meta$type=="static")
@@ -563,42 +569,81 @@ new_diffnet <- function(graph, toa, t0=min(toa, na.rm = TRUE), t1=max(toa, na.rm
 
 
   # Step 1.2: Checking that lengths fit
-  if (length(toa)!=meta$n) stop("-graph- and -toa- have different lengths (",
-                                meta$n, " and ", length(toa), " respectively). ",
-                                "-toa- should be of length n (number of vertices).")
-
-  # Step 2.1: Checking class of TOA and coercing if necesary -------------------
-  if (!inherits(toa, "integer")) {
-    warning("Coercing -toa- into integer.")
-    toa <- as.integer(toa)
+  if (num_of_behaviors == 1) {
+    if (length(toa)!=meta$n){ stop("-graph- and -toa- have different lengths (", meta$n, " and ", length(toa),
+         " respectively). ", "-toa- should be of length n (number of vertices).") }
+  } else {
+    if (length(toa[,1])!=meta$n) {stop("-graph- and -toa[,1]- have different lengths (", meta$n, " and ", length(toa[,1]),
+                                       " respectively). ", "-toa- should be of length n (number of vertices).") }
+  }
+
+  # Step 2.1: Checking class of TOA and coercing if necessary -------------------
+  if (num_of_behaviors==1) {
+    if (!inherits(toa, "integer")) {
+      warning("Coercing -toa- into integer.")
+      toa <- as.integer(toa)
+    }
+  } else {
+    for (q in 1:num_of_behaviors) {
+      if (!inherits(toa[,q], "integer")) {
+        warning("Coercing -toa- into integer.")
+        toa[,q] <- as.integer(toa[,q])
+      }
+    }
   }
 
   # Step 2.2: Checking names of toa
-  if (!length(names(toa)))
-    names(toa) <- meta$ids
+  if (num_of_behaviors==1) {
+    if (!length(names(toa))) {names(toa) <- meta$ids}
+  } else {
+    if (!length(rownames(toa))) { # Not necessary? toa_mat(toa, labels = meta$ids, t0=t0, t1=t1) already has labels
+      rownames(toa) <- meta$ids
+    }
+  }
 
   # Step 3.1: Creating Time of adoption matrix ---------------------------------
   mat <- toa_mat(toa, labels = meta$ids, t0=t0, t1=t1)
 
   # Step 3.2: Verifying dimensions and fixing meta$pers
 
-  if (meta$type != "static") {
-    tdiff <- meta$nper - ncol(mat[[1]])
-    if (tdiff < 0)
-      stop("Range of -toa- is bigger than the number of slices in -graph- (",
-           ncol(mat[[1]]), " and ", length(graph) ," respectively). ",
-           "There must be at least as many slices as range of toa.")
-    else if (tdiff > 0)
-      stop("Range of -toa- is smaller than the number of slices in -graph- (",
-           ncol(mat[[1]]), " and ", length(graph) ," respectively). ",
-           "Please provide lower and upper boundaries for the values in -toa- ",
-           "using -t0- and -t- (see ?toa_mat).")
+  if (num_of_behaviors==1) {
+    if (meta$type != "static") {
+      tdiff <- meta$nper - ncol(mat$adopt)
+      if (tdiff < 0)
+        stop("Range of -toa- is bigger than the number of slices in -graph- (",
+             ncol(mat$adopt), " and ", length(graph) ," respectively). ",
+             "There must be at least as many slices as range of toa.")
+      else if (tdiff > 0)
+        stop("Range of -toa- is smaller than the number of slices in -graph- (",
+             ncol(mat$adopt), " and ", length(graph) ," respectively). ",
+             "Please provide lower and upper boundaries for the values in -toa- ",
+             "using -t0- and -t- (see ?toa_mat).")
+    } else {
+      graph <- lapply(1:ncol(mat$adopt), function(x) methods::as(graph, "dgCMatrix"))
+      meta  <- classify_graph(graph)
+    }
   } else {
-    graph <- lapply(1:ncol(mat[[1]]), function(x) methods::as(graph, "dgCMatrix"))
-    meta  <- classify_graph(graph)
+    if (meta$type != "static") {
+      tdiff <- meta$nper - ncol(mat[[1]]$adopt)
+      if (tdiff < 0)
+        stop("Range of -toa- is bigger than the number of slices in -graph- (",
+             ncol(mat[[1]]$adopt), " and ", length(graph) ," respectively). ",
+             "There must be at least as many slices as range of toa.")
+      else if (tdiff > 0)
+        stop("Range of -toa- is smaller than the number of slices in -graph- (",
+             ncol(mat[[1]]$adopt), " and ", length(graph) ," respectively). ",
+             "Please provide lower and upper boundaries for the values in -toa- ",
+             "using -t0- and -t- (see ?toa_mat).")
+    } else {
+      graph <- lapply(1:ncol(mat[[1]]$adopt), function(x) methods::as(graph, "dgCMatrix"))
+      meta  <- classify_graph(graph)
+    }
   }
 
-  meta$pers <- as.integer(colnames(mat$adopt))
+  # labels of the time periods
+  if (num_of_behaviors==1) {
+    meta$pers <- as.integer(colnames(mat$adopt))
+  } else {meta$pers <- as.integer(colnames(mat[[1]]$adopt))} # same for all behaviors
 
   # Step 4.0: Checking the attributes ------------------------------------------
 
@@ -629,21 +674,40 @@ new_diffnet <- function(graph, toa, t0=min(toa, na.rm = TRUE), t1=max(toa, na.rm
   meta$multiple   <- multiple
   meta$name       <- ifelse(!length(name), "", ifelse(is.na(name), "",
                                                       as.character(name)))
-  meta$behavior   <- ifelse(!length(behavior), "", ifelse(is.na(behavior), "",
-                                                          as.character(behavior)))
   meta$version    <- utils::packageVersion("netdiffuseR")
 
   # Removing dimnames
   graph                  <- Map(function(x) Matrix::unname(x), x=graph)
-  dimnames(toa)          <- NULL
-  dimnames(mat$adopt)    <- NULL
-  dimnames(mat$cumadopt) <- NULL
+  #dimnames(toa)          <- NULL
+
+  if (num_of_behaviors==1) {
+    meta$behavior   <- ifelse(!length(behavior), "", ifelse(is.na(behavior), "",
+                                                            as.character(behavior)))
+    dimnames(mat$adopt)    <- NULL
+    dimnames(mat$cumadopt) <- NULL
+
+    adopt <- mat$adopt
+    cumadopt <- mat$cumadopt
+  } else {
+    meta$behavior <- paste(unlist(behavior), collapse = ", ")
+
+    for (q in 1:num_of_behaviors) {
+      dimnames(mat[[q]]$adopt)    <- NULL
+      dimnames(mat[[q]]$cumadopt) <- NULL
+    }
+    adopt <- list()
+    cumadopt <- list()
+    for (q in 1:num_of_behaviors) {
+      adopt[[q]] <- mat[[q]]$adopt
+      cumadopt[[q]] <- mat[[q]]$cumadopt
+    }
+  }
 
   return(structure(list(
     graph = graph,
     toa   = toa,
-    adopt = mat$adopt,
-    cumadopt = mat$cumadopt,
+    adopt = adopt,
+    cumadopt = cumadopt,
     # Attributes
     vertex.static.attrs = vertex.static.attrs,
     vertex.dyn.attrs    = vertex.dyn.attrs,
diff --git a/R/diffnet-methods.r b/R/diffnet-methods.r
index 109a0f9..95b420e 100644
--- a/R/diffnet-methods.r
+++ b/R/diffnet-methods.r
@@ -94,6 +94,22 @@ print.diffnet <- function(x, ...) {
                      paste(head(meta$ids, 8), collapse=", "),
                      ifelse(meta$n>8, ", ...", "") ,")")
 
+    # Computing prevalence for multi-diff
+    single <- !inherits(cumadopt, "list")
+    if (!single) {
+      prevalence_all <- character(length(cumadopt))
+      for (q in 1:length(cumadopt)) {
+        prevalence <- formatC(sum(cumadopt[[q]][,meta$nper]) / meta$n, digits = 2, format="f")
+        prevalence_all[q] <- prevalence
+      }
+      prevalence_all <- paste(prevalence_all, collapse = ", ")
+
+      num_of_behavior <- as.character(length(cumadopt))
+    } else {
+      prevalence <- formatC(sum(cumadopt[,meta$nper])/meta$n, digits = 2, format="f")
+      num_of_behavior <- "1"
+    }
+
     cat(
     "Dynamic network of class -diffnet-",
     paste(" Name               :", meta$name),
@@ -101,9 +117,12 @@ print.diffnet <- function(x, ...) {
     paste(" # of nodes         :", nodesl ),
     paste(" # of time periods  :", meta$nper, sprintf("(%d - %d)", meta$pers[1], meta$pers[meta$nper])),
     paste(" Type               :", ifelse(meta$undirected, "undirected", "directed")),
-    paste(" Final prevalence   :",
-          formatC(sum(cumadopt[,meta$nper])/meta$n, digits = 2, format="f")
-          ),
+    paste(" Num of behaviors   :", num_of_behavior),
+    if (single) {
+      paste(" Final prevalence   :", prevalence)
+    } else {
+      paste(" Prevalence         :", prevalence_all)
+    },
     paste(" Static attributes  :", vsa),
     paste(" Dynamic attributes :", vda),
     sep="\n"
@@ -194,44 +213,93 @@ summary.diffnet <- function(
     # x <-nelements/(meta$n * (meta$n-1))
   }))
 
-  # Computing moran's I
-  if (!skip.moran) {
+  # identify single-diff from multi-diff
+  single <- !inherits(object$cumadopt, "list")
 
-    m <- matrix(NA, nrow=length(slices), ncol=4,
-                dimnames = list(NULL, c("moran_obs", "moran_exp", "moran_sd", "moran_pval")))
+  # Computing moran's I
+  if (single) {
+    if (!skip.moran) {
+      m <- matrix(NA, nrow=length(slices), ncol=4,
+                  dimnames = list(NULL, c("moran_obs", "moran_exp", "moran_sd", "moran_pval")))
+
+      for (i in 1:length(slices)) {
+        # Computing distances
+        g <- approx_geodesic(object$graph[[slices[i]]], ...)
+        # Inverting it (only the diagonal may have 0)
+        g@x <- 1/g@x
+
+        m[i,] <- unlist(moran(object$cumadopt[,slices[i]], g))
+      }
 
-    for (i in 1:length(slices)) {
-      # Computing distances
-      g <- approx_geodesic(object$graph[[slices[i]]], ...)
+    # Computing new adopters, cumadopt and hazard rate
+    ad <- colSums(object$adopt[,slices,drop=FALSE])
+    ca <- t(cumulative_adopt_count(object$cumadopt))[slices,-3, drop=FALSE]
+    hr <- t(hazard_rate(object$cumadopt, no.plot = TRUE))[slices,,drop=FALSE]
+
+    # Left censoring
+    lc <- sum(object$toa == meta$pers[1], na.rm = TRUE)
+    rc <- sum(is.na(object$toa), na.rm=TRUE)
+
+    out <- data.frame(
+      adopt = ad,
+      cum_adopt = ca[,1],
+      cum_adopt_pcent = ca[,2],
+      hazard = hr,
+      density=d
+    )
 
-      # Inverting it (only the diagonal may have 0)
-      g@x <- 1/g@x
+    if (!skip.moran) {
+      out <- cbind(out, m)
+    }
 
-      m[i,] <- unlist(moran(object$cumadopt[,slices[i]], g))
+    if (no.print) return(out)
     }
-  }
-  # Computing adopters, cumadopt and hazard rate
-  ad <- colSums(object$adopt[,slices,drop=FALSE])
-  ca <- t(cumulative_adopt_count(object$cumadopt))[slices,-3, drop=FALSE]
-  hr <- t(hazard_rate(object$cumadopt, no.plot = TRUE))[slices,,drop=FALSE]
-
-  # Left censoring
-  lc <- sum(object$toa == meta$pers[1], na.rm = TRUE)
-  rc <- sum(is.na(object$toa), na.rm=TRUE)
-
-  out <- data.frame(
-    adopt = ad,
-    cum_adopt = ca[,1],
-    cum_adopt_pcent = ca[,2],
-    hazard = hr,
-    density=d
-  )
 
-  if (!skip.moran) {
-    out <- cbind(out, m)
-  }
+  } else {
+    if (!skip.moran) {
+      out_list <- list()
+      data_beh_list <- list()
+      for (q in 1:length(object$cumadopt)) {
+        m <- matrix(NA, nrow=length(slices), ncol=4,
+                    dimnames = list(NULL, c("moran_obs", "moran_exp", "moran_sd", "moran_pval")))
+
+        for (i in 1:length(slices)) {
+
+          g <- approx_geodesic(object$graph[[slices[i]]], ...)
+          g@x <- 1/g@x
+
+          m[i,] <- unlist(moran(object$cumadopt[[q]][,slices[i]], g))
+        }
+
+        # Computing new adopters, cumadopt and hazard rate
+        ad <- colSums(object$adopt[[q]][,slices,drop=FALSE])
+        ca <- t(cumulative_adopt_count(object$cumadopt[[q]]))[slices,-3, drop=FALSE]
+        hr <- t(hazard_rate(object$cumadopt[[q]], no.plot = TRUE))[slices,,drop=FALSE]
+
+        # Left censoring
+        lc <- sum(object$toa[,q] == meta$pers[1], na.rm = TRUE)
+        rc <- sum(is.na(object$toa[,q]), na.rm=TRUE)
 
-  if (no.print) return(out)
+        #data_beh_list[[q]] <- list(ad, ca, hr, lc, rc)
+
+        out <- data.frame(
+          adopt = ad,
+          cum_adopt = ca[,1],
+          cum_adopt_pcent = ca[,2],
+          hazard = hr,
+          density=d
+        )
+
+        if (!skip.moran) {
+          out <- cbind(out, m)
+        }
+
+        if (no.print) return(out)
+
+        out_list[[q]] <- out
+      }
+    }
+  }
 
   # Function to print data.frames differently
   header <- c(" Period "," Adopters "," Cum Adopt. (%) ",
@@ -253,31 +321,14 @@ summary.diffnet <- function(
       rule,"\n",sep="")
   cat(header,"\n")
   cat(hline, "\n")
-  for (i in 1:nrow(out)) {
-    cat(sprintf(
-      paste0("%",slen,"s", collapse=" "),
-      qf(meta$pers[slices[i]],0), qf(out[i,1],0),
-      sprintf("%s (%s)",
-        qf(out$cum_adopt[i],0),
-        qf(out$cum_adopt_pcent[i])
-        ),
-      ifelse(i==1, "-",qf(out$hazard[i])), qf(out$density[i]),
-      if (!skip.moran) {
-        if (is.nan(out$moran_sd[i]))
-          " - "
-        else
-          sprintf("%s (%s) %-3s",
-                  qf(out$moran_obs[i]),
-                  qf(out$moran_sd[i]),
-                  ifelse(out$moran_pval[i] <= .01, "***",
-                         ifelse(out$moran_pval[i] <= .05, "**",
-                                ifelse(out$moran_pval[i] <= .10, "*", ""
-                              )))
-                )
-      } else ""
-    ), "\n")
-  }
 
+  if (single) {
+    summary_diffnet_out_display(out, slen, meta, slices, qf, skip.moran)
+  } else {
+    for (q in 1:length(object$cumadopt)) {
+      summary_diffnet_out_display(out_list[[q]], slen, meta, slices, qf, skip.moran)
+    }
+  }
 
   # print(out, digits=2)
 
@@ -294,6 +345,32 @@ summary.diffnet <- function(
   invisible(out)
 }
 
+summary_diffnet_out_display <- function(out, slen, meta, slices, qf, skip.moran) {
+  for (i in 1:nrow(out)) {
+    cat(sprintf(
+      paste0("%",slen,"s", collapse=" "),
+      qf(meta$pers[slices[i]],0), qf(out[i,1],0),
+      sprintf("%s (%s)",
+              qf(out$cum_adopt[i],0),
+              qf(out$cum_adopt_pcent[i])
+      ),
+      ifelse(i==1, "-",qf(out$hazard[i])), qf(out$density[i]),
+
+      if (!skip.moran) {
+        if (is.nan(out$moran_sd[i])) {" - "}
+        else {sprintf("%s (%s) %-3s",
+                      qf(out$moran_obs[i]),
+                      qf(out$moran_sd[i]),
+                      ifelse(out$moran_pval[i] <= .01, "***",
+                             ifelse(out$moran_pval[i] <= .05, "**",
+                                    ifelse(out$moran_pval[i] <= .10, "*", ""
+                                    )))
+        )}
+      } else {""}
+    ), "\n")
+  }
+}
+
 #' Plot the diffusion process
 #'
 #' Creates a colored network plot showing the structure of the graph through time
diff --git a/R/rdiffnet.r b/R/rdiffnet.r
index 659075c..f007ef1 100644
--- a/R/rdiffnet.r
+++ b/R/rdiffnet.r
@@ -92,7 +92,7 @@
 #' }
 #'
 #' @examples
-#' # Asimple example -----------------------------------------------------------
+#' # A simple example -----------------------------------------------------------
 #' set.seed(123)
 #' z <- rdiffnet(100,10)
 #' z
@@ -111,45 +111,43 @@
 #' @name rdiffnet
 NULL
 
-rdiffnet_make_threshold <- function(x, n) {
+rdiffnet_make_threshold <- function(x, n, num_of_behaviors) {
 
-  # Using sapply to compute the threshold
-  if (inherits(x, "function")) {
-
-    thr <- sapply(1:n, x)
-
-  } else if ((length(x)==1) && is.numeric(x)) {
-
-    thr <- rep(x, n)
-
-  } else {
-    # Setting depending on class
-    if (any(class(x) %in% c("data.frame", "matrix"))) {
+  # Check if x is a matrix or array with correct dimensions
+  if (is.matrix(x) || is.array(x)) {
+    if (!all(dim(x) == c(n, num_of_behaviors))) {
+      stop("Incorrect threshold input in function -rdiffnet_make_threshold-. The matrix/array must have dimensions ", n, "x", num_of_behaviors, ".")
+    }
+    return(as.matrix(x)) # Return the matrix as-is
+  } else if (!is.list(x) && num_of_behaviors > 1) {
+    x <- rep(list(x), num_of_behaviors)
+  }
 
-      thr <- as.vector(as.matrix(x))
+  # Make a list, for single diffusion
+  if (!is.list(x) && num_of_behaviors==1) {
+    x <- list(x)
+  }
 
-      # Must match the length of n
-      if (length(thr) != n)
-        stop("Incorrect length for -threshold.dist- (",length(x),")",
-             ". It should be a vector of length ",n,".")
+  thr <- matrix(NA, nrow = n, ncol = num_of_behaviors)
 
-    } else if (is.vector(x)) {
+  for (q in seq_len(num_of_behaviors)) {
+    if (inherits(x[[q]], "function")) {
+      set.seed(123)
+      thr[, q] <- sapply(1:n, function(j) x[[q]]())
 
-      thr <- x
+    } else if (is.numeric(x[[q]]) && length(x[[q]]) == 1) {
 
-      # Must match the length of n
-      if (length(thr) != n)
-        stop("Incorrect length for -threshold.dist- (",length(x),")",
-             ". It should be a vector of length ",n,".")
+      thr[, q] <- rep(x[[q]], n)
 
-    } else {
+    } else if (is.vector(x[[q]]) && length(x[[q]]) == n) {
 
-      stop("-threshold.dist- must be either a numeric vector of length -n-, a numeric scalar, or a function.")
+      thr[, q] <- x[[q]]
 
+    } else if (is.vector(x[[q]]) && length(x[[q]]) != n) {
+      stop("Incorrect threshold input in function -rdiffnet_make_threshold-.")
     }
   }
-
-  thr
+  return(thr)
 }
 
 rdiffnet_check_seed_graph <- function(seed.graph, rgraph.args, t, n) {
@@ -316,14 +314,13 @@ rdiffnet <- function(
     seed.p.adopt   = 0.05,
     seed.graph     = "scale-free",
     rgraph.args    = list(),
-    rewire         = TRUE, #set TRUE originally
+    rewire         = TRUE,
     rewire.args    = list(),
     threshold.dist = runif(n),
     exposure.args  = list(),
     name           = "A diffusion network",
     behavior       = "Random contagion",
-    stop.no.diff   = TRUE,
-    behavior.num   = 1
+    stop.no.diff   = TRUE
   ) {
 
   # Checking options
@@ -335,6 +332,16 @@ rdiffnet <- function(
   if (!length(exposure.args[["valued"]])) exposure.args[["valued"]] <- getOption("diffnet.valued", FALSE)
   if (!length(exposure.args[["normalized"]])) exposure.args[["normalized"]] <- TRUE
 
+  if (inherits(exposure.args[["attrs"]], "matrix")) {
+    # Checking if the attrs matrix is has dims n x t
+    if (any(dim(exposure.args[["attrs"]]) != dim(matrix(NA, nrow = n, ncol = t)))) {
+      stop("Incorrect size for -attrs- in rdiffnet. Does not match n dim or t dim.")}
+    attrs_arr <- exposure.args[["attrs"]]
+    if (inherits(seed.p.adopt, "list")){
+      attrs_arr <- array(attrs_arr, dim = c(n, t, length(seed.p.adopt)))
+    } else {attrs_arr <- array(attrs_arr, dim = c(n, t, 1))}
+  }
+
   # Step 0.0: Creating the network seed ----------------------------------------
   # Checking the class of the seed.graph
   sgraph <- rdiffnet_check_seed_graph(seed.graph, rgraph.args, t, n)
@@ -376,44 +383,35 @@ rdiffnet <- function(
 
   # Step 1.0: Setting the seed nodes -----------------------------------------
 
-  # Step 1.1: Number of initial adopters
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
 
-  if (length(seed.p.adopt)>1 && length(seed.p.adopt) == behavior.num) {
+  seed.p.adopt <- rdiffnet_args$seed.p.adopt
+  seed.nodes <- rdiffnet_args$seed.nodes
+  behavior <- rdiffnet_args$behavior
+  num_of_behaviors <- rdiffnet_args$num_of_behaviors
 
-    n0 <- list()
-
-    for (i in seq_along(seed.p.adopt)) {
-
-      if ((seed.p.adopt[i] > 1) | (seed.p.adopt[i] < 0)) {
-        stop(paste("The proportion of initial adopters for behavior", i, "should be a number in [0,1]"))
-      }
-      if (n*seed.p.adopt[i] < 1) {
-        warning(paste("Set of initial adopters for behavior", i, "set to 1."))
-      }
+  # Step 1.1: Number of initial adopters
 
-      n0[[i]] <- max(1, n * seed.p.adopt[i])
-    }
+  n0 <- list()
 
-  } else if (length(seed.p.adopt)== 1 && behavior.num == 1) {
+  for (i in 1:num_of_behaviors) {
 
-    if ((seed.p.adopt > 1) | (seed.p.adopt < 0)) {
-      stop("The proportion of initial adopters should be a number in [0,1]")
+    if ((seed.p.adopt[[i]] > 1) | (seed.p.adopt[[i]] < 0)) {
+      stop(paste("The proportion of initial adopters for behavior", i, "should be a number in [0,1]"))
     }
-    if (n*seed.p.adopt < 1) {
-      warning("Set of initial adopters set to 1.")
+    if (n*seed.p.adopt[[i]] < 1) {
+      warning(paste("Set of initial adopters for behavior", i, "set to 1."))
     }
 
-    n0 <- max(1, n*seed.p.adopt)
-  } else {
-    stop("Error in setting number of initial adopters. Mismatch between length(seed.p.adopt) and behavior.num")
+    n0[[i]] <- max(1, n * seed.p.adopt[[i]])
   }
 
+  # Step 1.2: finding the nodes
+
+  d <- list()
 
-  # Step 1.2: Finding seed nodes
-  if (length(seed.nodes) > 1 && length(seed.nodes) == behavior.num && class(seed.nodes)!="list") {
-    # multi-diff. Something like seed.nodes <- c("marginal", "central"), and behavior.num <- 2
+  if (all(sapply(seed.nodes, is.character))) { # "central", "marginal", or "random"
 
-    d <- list()
     if (any(seed.nodes %in% c("central", "marginal"))) {
       dg <- dgr(sgraph)[, 1, drop = FALSE]
       central_d <- rownames(dg[order(dg, decreasing = TRUE), , drop = FALSE])
@@ -421,149 +419,227 @@ rdiffnet <- function(
     }
 
     for (i in seq_along(seed.nodes)) { # assign nodes characters values in seed.nodes
-      d[[i]] <- switch(seed.nodes[i],
-                      "central" = as.numeric(central_d[1:floor(n0[[i]])]),
-                      "marginal" = as.numeric(marginal_d[1:floor(n0[[i]])]),
-                      "random" = sample.int(n, floor(n0[[i]])),
-      stop("Unsupported -seed.nodes- value. It must be either \"central\", \"marginal\", or \"random\"")
-      )
+      d[[i]] <- switch(seed.nodes[[i]],
+                       "central" = as.numeric(central_d[1:floor(n0[[i]])]),
+                       "marginal" = as.numeric(marginal_d[1:floor(n0[[i]])]),
+                       "random" = sample.int(n, floor(n0[[i]])),
+                       stop("Unsupported -seed.nodes- value. It must be either \"central\", \"marginal\", or \"random\"")
+                      )
     }
-  } else if (length(seed.nodes) == 1 && behavior.num == 1) {
-    # Single-diff. Something like seed.nodes <- "central"
 
-    if (seed.nodes %in% c("central","marginal")) {
+  } else if (all(sapply(seed.nodes, is.numeric))) { # specific nodes
 
-      # Creating a degree ranking
-      d <- dgr(sgraph)[,1,drop=FALSE]
-      decre <- ifelse(seed.nodes == "central", TRUE, FALSE)
-      d <- rownames(d[order(d, decreasing = decre),,drop=FALSE])
-      d <- d[1:floor(n0)]
-      d <- as.numeric(d)
+    for (i in 1:num_of_behaviors) {
+      d[[i]] <- seed.nodes[[i]]
+    }
 
-    } else if (seed.nodes == "random") {
+  } else {
+    stop("Unsupported -seed.nodes- value. See the manual for references.")
+  }
 
-      d <- sample.int(n, floor(n0))
+  # Step 1.3: Defining cumadopt and toa (time of adoption) --------------------
 
-    } else {
-      stop("Unsupported -seed.nodes- value. It must be either \"central\", \"marginal\", or \"random\"")
+  cumadopt <- array(0L, dim = c(n, t, num_of_behaviors))
+
+  toa <- matrix(NA, nrow = dim(cumadopt)[1], ncol = dim(cumadopt)[3])
+
+  for (i in 1:num_of_behaviors) {
+    cumadopt[d[[i]],,i] <- 1L
+  }
+
+  # Step 2.0: Thresholds -------------------------------------------------------
+
+  thr <- rdiffnet_make_threshold(threshold.dist, n, num_of_behaviors)
+
+  # Step 3.0: Running the simulation -------------------------------------------
+
+  for (i in 2:t) {
+    if (exists("attrs_arr")){
+      exposure.args[c("attrs")] <- list(attrs_arr[,i, ,drop=FALSE])
     }
-  } else if (class(seed.nodes)=="list" && length(seed.nodes) != behavior.num) {
-    # Something like seed.nodes <- c("marginal", "central"), BUT behavior.num <- 3
-    stop("Error in finding seed nodes. Mismatch between length(seed.nodes) and behavior.num")
 
-  } else if (!inherits(seed.nodes, "character")) {
+    exposure.args[c("graph", "cumadopt")] <- list(sgraph[i], cumadopt[,i, ,drop=FALSE])
+    expo <- do.call(exposure, exposure.args)
+
+    for (q in 1:num_of_behaviors) {
+
+      whoadopts <- which( (expo[,,q] >= thr[,q]) )
+      cumadopt[whoadopts, i:t, q] <- 1L
+                                                        # ADD SOMETHING TO DISADOPT
+
+      toa[, q] <- apply(cumadopt[,, q], 1, function(x) {
+        first_adopt <- which(x == 1)
+        if (length(first_adopt) > 0) first_adopt[1] else NA
+      })
 
-    if (class(seed.nodes)=="list" && length(seed.nodes) != behavior.num) {
-      # Something like seed.nodes <- list(c(1,4), c(3,6,8)), BUT behavior.num <- 3
-      stop("Particular seed nodes provided. Mismatch between length(seed.nodes) and behavior.num")
-    } else {
-      # single-diff and multi-diff. # Something like seed.nodes <- c(3,6,8)), AND behavior.num <- 1,
-      # or seed.nodes <- list(c(1,4), c(3,6,8)), AND behavior.num <- 2
-      d <- seed.nodes
     }
-  } else {stop("Unsupported -seed.nodes- value. See the manual for references.") }
+  }
 
-  # Step 1.3: Defining cumadopt and toa (time of adoption) --------------------
+  for (i in 1:num_of_behaviors) {
+    reachedt <- max(toa[,i], na.rm=TRUE)
+
+    if (reachedt == 1) {
+      if (stop.no.diff)
+        stop(paste("No diffusion in this network for behavior", i, "(Ups!) try changing the seed or the parameters."))
+      else
+        warning(paste("No diffusion for behavior", i, " in this network."))
+    }
+  }
 
-  if (class(d) == "list") {
-    # multi-diff
+  # Step 4.0: Creating diffnet object ------------------------------------------
+  # Checking attributes
+  isself <- any(sapply(sgraph, function(x) any(Matrix::diag(x) != 0) ))
 
-    if (length(d) != behavior.num) {
-      stop("Error: length(d) must be the same as behavior.num")
+  if (num_of_behaviors==1) {
+    toa <- as.integer(toa)
     }
 
-    cumadopt <- array(0L, dim = c(n, t, behavior.num))
+  new_diffnet(
+    graph      = sgraph,
+    toa        = toa,
+    self       = isself,
+    t0         = 1,
+    t1         = t,
+    vertex.static.attrs = data.frame(real_threshold=thr),
+    name       = name,
+    behavior   = behavior
+  )
+}
+
+
+rdiffnet_validate_args <- function(seed.p.adopt, seed.nodes, behavior) {
+
+  # seed.p.adopt stuff
+
+  # The class of seed.p.adopt determines if is a single or multiple diff pross.
+
+  if (inherits(seed.p.adopt, "list")) {
+
+    message(paste("Message: Multi-diffusion behavior simulation selected.",
+                  "Number of behaviors: ", length(seed.p.adopt)))
+
+    multi <- TRUE
 
-    # Setting seed nodes via array
-    for (i in seq_along(d)) {
-      cumadopt[d[[i]],,i] <- 1L
+  } else if (inherits(seed.p.adopt, "numeric")) {
+
+    if (length(seed.p.adopt)>1) {
+      stop(paste("length(seed.p.adopt) =", length(seed.p.adopt),
+                 ", but for multi-diffusion -seed.p.adopt- must be a -list-."))
     }
+
+    multi <- FALSE
+
   } else {
-    # single-diff
-    cumadopt <- matrix(0L, ncol=t, nrow=n)
-    toa      <- matrix(NA, ncol=1, nrow= n)
 
-    # Setting seed nodes via vector
-    toa[d]       <- 1L                # REMINDER TO DELETE THIS OBJECT !!!
-    cumadopt[d,] <- 1L
+    stop("The object -seed.p.adopt- must be a -numeric- (for a single behavior diff)",
+         "or a -list- (multiple behavior diff).")
   }
 
-  # Step 2.0: Thresholds -------------------------------------------------------
-  thr <- rdiffnet_make_threshold(threshold.dist, n) # REMINDER TO CHANGE rdiffnet_make_threshold
+  # seed.nodes stuff
 
-  # Step 3.0: Running the simulation -------------------------------------------
-  for (i in 2:t) {
-    if (!is.na(dim(cumadopt)[3])) {
-      # multi-diff. Computing exposure
-                                                    # ONLY MEANWHILE
-      thr <- array(c(thr,rev(thr)), dim=c(length(thr), dim(cumadopt)[3]))
-
-      exposure.args[c("graph", "cumadopt")] <- list(sgraph[i], cumadopt[,i, ,drop=FALSE])
-      expo <- do.call(exposure, exposure.args)
-      #for (q in 1:dim(cumadopt)[3]) {
-      #  exposure.args[c("graph", "cumadopt")] <- list(sgraph[i], cumadopt[,i,q,drop=FALSE])
-      #}
-
-      toa <- matrix(NA, nrow = dim(cumadopt)[1], ncol = dim(cumadopt)[3])
-
-      for (q in 1:dim(cumadopt)[3]) {
-        whoadopts <- which( (expo[,,q] >= thr[,q]) )# & is.na(toa))
-        cumadopt[whoadopts, i:t, q] <- 1L
-                                                      # ADD SOMETHING TO DISADOPT
-        # Initialize 'toa' with NA values
-        toa[, q] <- apply(cumadopt[,, q], 1, function(x) {
-          first_adopt <- which(x == 1)
-          if (length(first_adopt) > 0) first_adopt[1] else NA
-        })
+  if (multi) {
+
+    # For multi-diff.
+
+    if (inherits(seed.nodes, "list")) {
+      if (length(seed.nodes) != length(seed.p.adopt)) {
+        stop("Length of lists -seed.nodes- and -seed.p.adopt- must be the same for multi diffusion.")
       }
 
-    } else {
-      # single-diff. Computing exposure
-      exposure.args[c("graph", "cumadopt")] <- list(sgraph[i], cumadopt[,i,drop=FALSE])
-      expo <- do.call(exposure, exposure.args)
+      if (all(sapply(seed.nodes, is.character))) {
+
+        if (any(!seed.nodes %in% c("marginal", "central", "random"))) {
+          stop("Some element in list -seed.nodes- is a -character- different from 'marginal', 'central', or 'random'.")
+        }
+
+      } else if (all(sapply(seed.nodes, is.numeric))) {
+
+        if (any(sapply(seed.nodes, is.null))) {
+          stop("There is a NULL -numeric- element")
+        }
 
-      whoadopts <- which( (expo >= thr) & is.na(toa))
-      toa[whoadopts] <- i
-      cumadopt[whoadopts, i:t] <- 1L
+        if (any(sapply(seed.nodes, function(x) any(x != round(x))))) {
+          stop("Some value in the elements of the list -seed.nodes- is non-integer.")
+        }
+
+      } else {
+        stop("All elements of the list seed.nodes must be either -character- or -numeric-.")
+      }
+    } else if (inherits(seed.nodes, "numeric")) {
+      message("Message: Object -seed.nodes- converted to a -list-.",
+              "All behaviors will have the same -", seed.nodes, "- seed nodes.")
+
+      seed.nodes <- replicate(length(seed.p.adopt), seed.nodes, simplify = FALSE)
+    } else if (inherits(seed.nodes, "character")) {
+      if (length(seed.nodes)==length(seed.p.adopt)) {
+        seed.nodes <- as.list(seed.nodes)
+        message("Message: Object -seed.nodes- converted to a -list-.",
+                "For example, the first behavior has seed -", seed.nodes[[1]], "-, the second has -", seed.nodes[[2]], "-, etc.")
+      } else {
+
+      message("Message: Object -seed.nodes- converted to a -list-.",
+              "All behaviors will have the same -", seed.nodes, "- seed nodes.")
+      seed.nodes <- replicate(length(seed.p.adopt), seed.nodes, simplify = FALSE)
+      }
+    } else {
+      stop("Unsupported -seed.nodes- value. See the manual for references.")
     }
-  }
-                                      # GENERALIZE TO MULTI-DIFF
-  reachedt <- max(toa, na.rm=TRUE)
 
-  # Checking the result
-  if (reachedt == 1) {
-    if (stop.no.diff)
-      stop("No diffusion in this network (Ups!) try changing the seed or the parameters.")
-    else
-      warning("No diffusion in this network.")
-  }
+    if (inherits(behavior, "list")) {
+      if (length(seed.p.adopt)!=length(behavior)) {
+        stop("If -behavior- is a list, it must be of the same length as -seed.p.adopt-.")
+      }
+    } else if (inherits(behavior, "character") && length(behavior) > 1) {
+      if (length(behavior) != length(seed.p.adopt)) {
+        stop("Mismatch between length(behavior) and length(seed.p.adopt)")
+      } else {
+        behavior <- as.list(behavior)
+      }
+    } else if (inherits(behavior, "character") && length(behavior) == 1) {
+      message(paste("Message: Name of 1 behavior provided, but", length(seed.p.adopt), "are needed. "),
+              "Names generalized to 'behavior'_1, 'behavior'_2, etc.")
+      behaviors <- list()
 
-  # Step 4.0: Creating diffnet object ------------------------------------------
-  # Checking attributes
-  isself <- any(sapply(sgraph, function(x) any(Matrix::diag(x) != 0) ))
+      for (i in seq_along(seed.p.adopt)) {
+        behaviors[[i]] <- paste(behavior, i, sep = "_")
+      }
+
+      behavior <- behaviors
+    }
 
-  if (!is.na(dim(cumadopt)[3])) {
-    new_diffnet(
-      graph      = sgraph,
-      toa        = toa,
-      self       = isself,
-      t0         = 1,
-      t1         = t,
-      vertex.static.attrs = data.frame(real_threshold=thr),
-      name       = name,
-      behavior   = behavior
-    )
   } else {
-    new_diffnet(
-      graph      = sgraph,
-      toa        = as.integer(toa),
-      self       = isself,
-      t0         = 1,
-      t1         = t,
-      vertex.static.attrs = data.frame(real_threshold=thr),
-      name       = name,
-      behavior   = behavior
-    )
+
+    # For Single-diff.
+
+    if (length(seed.nodes) == 1 && inherits(seed.nodes, "character")) {
+
+      if (!seed.nodes %in% c("marginal", "central", "random")) {
+        stop("Object -seed.nodes- is a -character- different from 'marginal', 'central', or 'random'.")
+      }
+
+    } else if (!inherits(seed.nodes, "character")) {
+
+      if (any(sapply(seed.nodes, function(x) any(x != round(x))))) {
+        stop("Some value in the elements of the list -seed.nodes- is non-integer.")
+      }
+
+    } else {
+      stop("Unsupported -seed.nodes- value. See the manual for references.")
+    }
+
+    if (length(behavior)>1) {
+      stop("More names were provided than necessary.")
+    }
+
+    seed.p.adopt <- list(seed.p.adopt)
+    seed.nodes <- list(seed.nodes)
+    behavior <- list(behavior)
   }
-}
 
+  list(
+    seed.p.adopt = seed.p.adopt,
+    seed.nodes = seed.nodes,
+    behavior = behavior,
+    num_of_behaviors = length(seed.p.adopt)
+  )
+}
diff --git a/R/rewire.r b/R/rewire.r
index d46479d..2cecbd3 100644
--- a/R/rewire.r
+++ b/R/rewire.r
@@ -234,7 +234,7 @@ rewire_graph <- function(
 
   if (copy.first) {
 
-    warning(
+    message(
       "The option -copy.first- is set to TRUE. In this case, the first graph will be ",
       "treated as a baseline, and thus, networks after T=1 will be replaced with T-1.",
       immediate. = TRUE
@@ -407,7 +407,7 @@ rewire_graph.array <-function(
 #'
 #' Mantel, N. (1967). The detection of disease clustering and a generalized
 #' regression approach. Cancer Research, 27(2), 209–20.
-#' 
+#'
 #' @seealso This function can be used as null distribution in \code{struct_test}
 #' @family simulation functions
 #' @export
diff --git a/R/stats.R b/R/stats.R
index 501f10a..8149b2f 100644
--- a/R/stats.R
+++ b/R/stats.R
@@ -485,7 +485,7 @@ NULL
       if (normalized) {
         ans[,q] <- as.vector(graph %*% (attrs * cumadopt[,,q]) / norm)
       } else {
-        ans[,q] <- as.vector(graph_slice %*% (attrs * cumadopt[,,q]))
+        ans[,q] <- as.vector(graph %*% (attrs * cumadopt[,,q]))
       }
     }
   } else {
@@ -566,7 +566,11 @@ exposure <- function(
 
   # Checking attrs
   if (!length(attrs)) {
-    attrs <- matrix(1, ncol=ncol(cumadopt), nrow=nrow(cumadopt))
+    if (!is.na(dim(cumadopt)[3])) {
+    attrs <- array(1, dim = c(nrow(cumadopt), ncol(cumadopt), 1))}
+    else {attrs <- matrix(1, ncol=ncol(cumadopt), nrow=nrow(cumadopt))}
+  } else if (!is.na(dim(cumadopt)[3])) {
+    attrs <- array(attrs, dim = c(nrow(attrs), ncol(attrs), 1))
   }
 
   # Checking alt graph
@@ -626,20 +630,18 @@ exposure.list <- function(
   #  degree, indegree, outdegree, or a user defined vector.
   #  by default is user equal to 1
 
-  dim_attrs <- dim(attrs) # default n x T matrix of 1's
-  if (!length(dim_attrs)) stop("-attrs- must be a matrix of size n by T.")
+  # dim(attrs) default n x T matrix of 1's
+  if (!length(dim(attrs))) stop("-attrs- must be a matrix of size n by T.")
 
   if (!is.na(dim(cumadopt)[3])) {
-    attrs_mul <- array(rep(attrs, dim(cumadopt)[3]), dim = c(dim_attrs, dim(cumadopt)[3]))
-    dim_attrs <- dim(attrs_mul) # now n x T x q array of 1's, q behaviors
-    if (any(dim_attrs != dim(cumadopt))) stop("Incorrect size for -attrs-. ",
+    if (dim(cumadopt)[3]>1 && any(dim(attrs)[-3] != dim(cumadopt)[-3])) stop("Incorrect size for -attrs-. ",
                                               "Does not match n dim or t dim.")
   } else {
-    if (any(dim_attrs != dim(cumadopt))) stop("Incorrect size for -attrs-. ",
+    if (any(dim(attrs) != dim(cumadopt))) stop("Incorrect size for -attrs-. ",
                                               "It must be of size that -cumadopt-.")
   }
-
   add_dimnames.mat(cumadopt)
+  add_dimnames.mat(attrs)
 
   output <- exposure_for(graph, cumadopt, attrs, outgoing, valued, normalized,
                          self, lags)
@@ -667,7 +669,7 @@ exposure_for <- function(
       for (i in 1:(nslices(graph) - lags)) {
         out[, i + lags, ] <- .exposure(graph[[i]],
                                        cumadopt[, i, , drop = FALSE],
-                                       attrs[, i, drop = FALSE],
+                                       attrs[, i, , drop = FALSE],
                                        outgoing = outgoing,
                                        valued = valued,
                                        normalized = normalized,
@@ -677,7 +679,7 @@ exposure_for <- function(
       for (i in (1 - lags):nslices(graph)) {
         out[, i + lags, ] <- .exposure(graph[[i]],
                                        cumadopt[, i, , drop = FALSE],
-                                       attrs[, i, drop = FALSE],
+                                       attrs[, i, , drop = FALSE],
                                        outgoing = outgoing,
                                        valued = valued,
                                        normalized = normalized,
diff --git a/tests/testthat/test-adjmat.R b/tests/testthat/test-adjmat.R
index a037278..1cf02bb 100644
--- a/tests/testthat/test-adjmat.R
+++ b/tests/testthat/test-adjmat.R
@@ -145,6 +145,69 @@ test_that("Checking toa_mat output", {
   expect_equal(amat, toa_mat(diffnet))
 })
 
+################################################################################
+# Time of adoption (multiple)
+################################################################################
+context("Time of Adoption -multiple- (toa_mat, toa_dif)")
+
+times_1 <- c(2001L, 2004L, 2003L, 2008L)
+times_2 <- c(2001L, 2005L, 2006L, 2008L)
+times <- matrix(c(times_1, times_2), nrow = 4, ncol = 2)
+graph <- lapply(2001:2008, function(x) rgraph_er(4))
+diffnet <- new_diffnet(graph, times)
+#toa_mat(diffnet)
+toa <- toa_mat(times)
+
+test_that("Should warn about -times- not been integer. -multiple-.", {
+  times_1 <- as.numeric(c(2001L, 2004L, 2003L, 2008L))
+  times_2 <- as.numeric(c(2001L, 2002L, 2003L, 2005L))
+  times <- matrix(c(times_1, times_2), nrow = 4, ncol = 2)
+  expect_warning(toa_mat(times), "will be coersed to integer")
+})
+
+test_that("Dimensions of TOA mat should be ok. -multiple-.", {
+  for (q in 1:length(toa)) {
+    cumadopt <- t(apply(toa[[q]]$adopt, 1, cumsum))
+    expect_equal(dim(toa[[q]]$adopt), c(4, length(min(times[,q]):max(times[,q]))))
+    expect_equal(dim(toa[[q]]$adopt), dim(toa[[q]]$cumadopt), info = "adopt and cumadopt are equal dim")
+    expect_equal(t(apply(toa[[q]]$adopt, 1, cumsum)), toa[[q]]$cumadopt, info = "cumadopt is the cumsum")
+  }
+})
+
+test_that("Passing labels should work. -multiple-.", {
+  for (q in 1:length(toa)) {
+    labs <- letters[1:length(times[,q])]
+    toa_q <- toa_mat(times[,q], labels=labs)
+    expect_equal(rownames(toa_q$adopt), labs)
+    expect_equal(rownames(toa_q$cumadopt), labs)
+  }
+})
+
+# test_that("In toa_diff, its dim should be equal to the input mat. -multiple-.", {
+#   expect_equal(dim(toa_diff(times)), c(4,4))
+#   expect_equal(dim(toa_diff(as.integer(times))), c(4,4))
+#   expect_equal(toa_diff(times), toa_diff(as.integer(times)))
+# })
+
+test_that("Checking toa_mat output. -multiple-.", {
+
+  # Manual calc
+  mat <- matrix(0, nrow=4, ncol=8)
+  mat <- array(rep(mat,2), dim = c(nrow(mat), ncol(mat), 2))
+  dimnames(mat) <- list(1:4, 2001:2008)
+  amat_tot <- list()
+  for (q in 1:dim(mat)[3]) {
+    amat <- list(adopt=mat[,,q], cumadopt=mat[,,q])
+    for (i in 1:4) {
+      amat$adopt[i,times[i,q] - 2000] <- 1
+      amat$cumadopt[i,] <- cumsum(amat$adopt[i,])
+    }
+    amat_tot[[q]] <- amat
+  }
+
+  expect_equal(amat_tot, toa_mat(diffnet))
+})
+
 ################################################################################
 # Isolated
 ################################################################################
diff --git a/tests/testthat/test-rdiffnet-parameters.R b/tests/testthat/test-rdiffnet-parameters.R
new file mode 100644
index 0000000..7fa3e18
--- /dev/null
+++ b/tests/testthat/test-rdiffnet-parameters.R
@@ -0,0 +1,224 @@
+# Single --------------------------------------------------------------------
+
+# Must work
+test_that(
+  "Checking single diffusion rdiffnet args", {
+
+  # Must work
+
+  seed.p.adopt <- c(0.14)
+  seed.nodes <- c('random')
+  behavior <- c("random behavior")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+
+  expect_type(rdiffnet_args$seed.p.adopt, "list")
+  expect_type(rdiffnet_args$seed.nodes, "list")
+  expect_type(rdiffnet_args$behavior, "list")
+
+  seed.p.adopt <- 0.14
+  seed.nodes <- 'random'
+  behavior <- "random behavior"
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+
+  expect_type(rdiffnet_args$seed.p.adopt, "list")
+  expect_type(rdiffnet_args$seed.nodes, "list")
+  expect_type(rdiffnet_args$behavior, "list")
+
+  # Must show ERROR
+
+  seed.p.adopt <- c(0.4,0.82)
+  seed.nodes <- c('random')
+  behavior <- "random behavior"
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+
+  seed.p.adopt <- c(0.14)
+  seed.nodes <- c('random', 'central')
+  behavior <- "random behavior"
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+
+  seed.p.adopt <- c(0.14)
+  seed.nodes <- "central"
+  behavior <- c("random behavior_1", "random behavior_2")
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+  behavior <- list("random behavior_1", "random behavior_2")
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+})
+
+test_that("Checking threshold for single diffusion", {
+
+  n <- 50
+  num_of_behaviors <- 1
+
+  # Must work
+
+  x <- 0.35 # numeric scalar
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n, 1))
+
+  x <- runif(n) # vector of length n
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n, 1))
+
+  x <- function() runif(1) # function
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n, 1))
+
+
+  # Must show ERROR
+
+  x <- runif(100)# Length greater than n
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors),
+    "Incorrect threshold input in function -rdiffnet_make_threshold-."
+  )
+
+  x <- runif(25)# Length less than n
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors),
+    "Incorrect threshold input in function -rdiffnet_make_threshold-."
+  )
+
+  x <- "invalid_input"# Non-numeric input
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors),
+    "Incorrect threshold input in function -rdiffnet_make_threshold-."
+  )
+
+})
+
+# Multiple --------------------------------------------------------------------
+
+test_that("Multi diff models rdiff args work", {
+
+  # Must work
+
+  seed.p.adopt <- list(0.14,0.05)
+  seed.nodes <- "random"
+  behavior <- "random behavior"
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$seed.p.adopt, "list")
+  expect_type(rdiffnet_args$seed.nodes, "list")
+  expect_type(rdiffnet_args$behavior, "list")
+
+  seed.nodes <- c(1,3,5)
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$seed.nodes, "list")
+
+  seed.nodes <- c('marginal',"central")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$seed.nodes, "list")
+
+  behavior <- list("random behavior_1", "random behavior_2")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$behavior, "list")
+
+  behavior <- c("random behavior_1", "random behavior_2")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$behavior, "list")
+
+  behavior <- "random behavior" #Default
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$behavior, "list")
+
+  behavior <- c("random behavior_1")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$behavior, "list")
+
+  seed.nodes <- c(1,3,5)
+  behavior <- list("random behavior_1", "random behavior_2")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$seed.nodes, "list")
+
+  seed.nodes <- list('marginal',"central")
+  rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  expect_type(rdiffnet_args$seed.nodes, "list")
+
+  # Must show ERROR
+
+  seed.p.adopt <- c(0.14,0.05)
+  seed.nodes <- list('random', "central")
+  behavior <- list("random behavior_1", "random behavior_2")
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+
+  behavior <- list("random behavior_1")
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+
+  seed.nodes <- list('random')
+  behavior <- list("random behavior_1", "random behavior_2")
+  expect_error(
+    rdiffnet_args <- rdiffnet_validate_args(seed.p.adopt, seed.nodes, behavior)
+  )
+})
+
+test_that("Checking threshold for multiple diffusion", {
+
+  n <- 50
+  num_of_behaviors <- 2
+
+  # Must work
+
+  # not list entries
+
+  x <- 0.35 # numeric scalar
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equivalent(thr, matrix(x, nrow=n, ncol=num_of_behaviors))
+
+  x <- runif(n) # vector of length n
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equivalent(thr, matrix(rep(x, num_of_behaviors), nrow = n, ncol = num_of_behaviors))
+
+  x <- function() runif(1) # function
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  set.seed(123)
+  expect_equal(thr, t(sapply(1:n, function(i) rep(x(), num_of_behaviors))))
+
+  # list entries
+
+  x <- matrix(runif(100), nrow = n, ncol = num_of_behaviors) # matrix input
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n, num_of_behaviors))
+
+  x <- list(function() runif(1), function() rexp(1)) # list of functions
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n, num_of_behaviors))
+
+  x <- list(0.14,0.05) # list of scalars
+  thr <- rdiffnet_make_threshold(x, n = n, num_of_behaviors = num_of_behaviors)
+  expect_equal(dim(thr), c(n,num_of_behaviors))
+
+  x <- list(runif(n), runif(n)) # list of vectors
+  thr <- rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors)
+  expect_equal(dim(thr),c(n,num_of_behaviors))
+
+
+  # Must show ERROR
+
+  x <- list(runif(2*n),runif(n)) # incorrect vector length (too long for one behavior)
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors),
+    "Incorrect threshold input in function -rdiffnet_make_threshold-."
+  )
+
+  x <- list(runif(n/2),runif(n)) # incorrect vector length (too short for one behavior)
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors),
+    "Incorrect threshold input in function -rdiffnet_make_threshold-."
+  )
+
+  x <- c(runif(n),runif(n)) # the input should be a list
+  expect_error(
+    rdiffnet_make_threshold(x,n=n,num_of_behaviors=num_of_behaviors)
+  )
+})
diff --git a/tests/testthat/test-rdiffnet.R b/tests/testthat/test-rdiffnet.R
index fb703e0..380ffeb 100644
--- a/tests/testthat/test-rdiffnet.R
+++ b/tests/testthat/test-rdiffnet.R
@@ -88,7 +88,7 @@ test_that("Error and warning on rdiffnet", {
 
   set.seed(111)
 
-  expect_error(rdiffnet(100, 5, threshold.dist = rep(10,10)), "Incorrect length")
+  expect_error(rdiffnet(100, 5, threshold.dist = rep(10,10)))
   expect_error(rdiffnet(100, 5, threshold.dist = rep(10,100)), "No diffusion")
   expect_warning(rdiffnet(100, 5, threshold.dist = rep(10,100), stop.no.diff = FALSE), "No diffusion")
 
@@ -106,3 +106,60 @@ test_that("Simulation study", {
   expect_equal(ans0, ans1)
 
 })
+
+# Testing diffnet class across several inputs (single)
+test_that("rdiffnet must run across several inputs (single)", {
+  expect_s3_class(rdiffnet(100, 5), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = 0.1), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = 0.1, seed.nodes = 'random'), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.nodes = c(1, 3, 5)), "diffnet")
+
+  # summary
+  net_1 <- rdiffnet(100, 5, seed.nodes = c(1,3,5))
+  expect_s3_class(summary(net_1), "data.frame")
+})
+
+# Testing diffnet class across several inputs (multiple)
+test_that("rdiffnet must run across several inputs (multiple)", {
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08)), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), behavior = c('tabacco', 'alcohol')), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), seed.nodes = 'random'), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), seed.nodes = c('random', 'central')), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = 0.3), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = list(0.1, 0.2)), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = rexp(100)), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = list(rexp(100), runif(100))), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = function(x) 0.3), "diffnet")
+  expect_s3_class(rdiffnet(100, 5, seed.p.adopt = list(0.1, 0.08), threshold.dist = list(function(x) 0.3, function(x) 0.2)), "diffnet")
+
+  net_2 <- rdiffnet(100, 5, seed.p.adopt = list(0.05,0.05), seed.nodes = c(1,3,5))
+  expect_s3_class(summary(net_2), "data.frame")
+})
+
+test_that("All should be equal! (multiple)", {
+  set.seed(12131)
+  n            <- 50
+  t            <- 5
+  graph        <- rgraph_ws(n, 4, p=.3)
+  seed.p.adopt <- list(0.1, 0.1)
+  seed.nodes   <- c(1,5,7,10)
+  thr          <- runif(n, .2,.4)
+  thr_list     <- list(thr,thr)
+
+  # Generating identical networks
+  net1 <- rdiffnet(seed.graph = graph, seed.nodes = seed.nodes, seed.p.adopt = seed.p.adopt,
+                   t = t, rewire = FALSE, threshold.dist = thr_list)
+
+  net2 <- rdiffnet(seed.graph = graph, seed.nodes = seed.nodes, seed.p.adopt = seed.p.adopt,
+                   t = t, rewire = FALSE, threshold.dist = thr_list)
+
+  expect_equal(net1, net2)
+})
+
+
+#rdiffnet(100, 5, seed.p.adopt = 0.9, threshold.dist = 2, exposure.args = list(normalized=FALSE))
+
+# set.seed(1234)
+# net1 <- rdiffnet(100, 5, rewire = FALSE, seed.p.adopt = list(0.1,0.08), seed.nodes = c(1,3,5))
+# net2 <- rdiffnet(100, 5, rewire = FALSE, seed.p.adopt = list(0.1,0.08), seed.nodes = list(c(1,3,5),c(1,3,5)))
+# expect_equal(net1, net2)