diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..2c46dca
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,11 @@
+^\.Rproj$
+^\.Rproj\.user$
+^Makefile$
+^CLAUDE\.md$
+^AGENTS\.md$
+^\.claude$
+^\.semquery$
+^\.beads$
+^local_data$
+^\.git$
+^\.gitignore$
diff --git a/.gitignore b/.gitignore
index af8a293..c4fbb93 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,4 @@
 *.o
 *.so
 *.dll
+tests/*.pdf
diff --git a/DESCRIPTION b/DESCRIPTION
index 176c04d..8e343f2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,6 +22,7 @@ RoxygenNote: 7.3.3
 Depends: R (>= 4.3.0)
 Imports:
     igraph,
+    magrittr,
     plotly,
     visNetwork,
     matrixStats,
@@ -33,7 +34,9 @@ Suggests:
     testthat (>= 3.0.0),
     knitr,
     rmarkdown,
-    BiocStyle
+    BiocStyle,
+    Rtsne,
+    uwot
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Config/roxygen2/version: 8.0.0
diff --git a/Makefile b/Makefile
index c53b20b..f80e9cd 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-PKGNAME := WGCNAplus
+PKGNAME := lasagna
 
 .PHONY: doc build install test check clean
 
diff --git a/NAMESPACE b/NAMESPACE
index 37447c9..076341f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,6 +14,7 @@ importFrom(grDevices,adjustcolor)
 importFrom(grDevices,rainbow)
 importFrom(graphics,plot)
 importFrom(graphics,text)
+importFrom(magrittr,"%>%")
 importFrom(matrixStats,rowSds)
 importFrom(plotly,add_markers)
 importFrom(plotly,add_mesh)
diff --git a/R/graph.R b/R/graph.R
index 066ed78..43abf68 100644
--- a/R/graph.R
+++ b/R/graph.R
@@ -43,61 +43,127 @@ prune_graph <- function(graph,
                         select = NULL,
                         prune = TRUE) {
 
-  if (is.null(layers)) layers <- graph$layers
-  if (is.null(layers)) layers <- unique(igraph::V(graph)$layer)
-  layers <- setdiff(layers, c("SOURCE", "SINK"))
+  if (is.null(graph)) {
+    message("[lasagna::prune_graph]: graph is NULL. Please check run of create_model() and solve(). Exiting.")
+    return(NULL)
+  }
+
+  layers <- resolve_layers(graph, layers)
   graph <- igraph::subgraph(graph, igraph::V(graph)$layer %in% layers)
 
   if (!"value" %in% names(igraph::vertex_attr(graph))) {
     stop("vertex must have 'value' attribute")
   }
 
-  ## select nodes/modules
-  if (!is.null(select)) {
-    v1 <- (igraph::V(graph)$name %in% select)
-    v2 <- (sub(".*:", "", igraph::V(graph)$name) %in% select)
-    v3 <- (sub(":.*", "", igraph::V(graph)$name) %in% select)
-    graph <- igraph::subgraph(graph, vids = which(v1 | v2 | v3))
-  }
-
   if (!is.null(filter)) {
     if (!is.list(filter)) stop("filter must be a named list")
     if (is.null(names(filter))) stop("filter must be a named list")
-    for (k in names(filter)) {
-      vv <- igraph::V(graph)$name
-      filt <- filter[[k]]
-      vids <- igraph::V(graph)$layer != k | grepl(filt, vv, ignore.case = TRUE)
-      graph <- igraph::subgraph(graph, which(vids))
-    }
   }
 
-  ## select ntop features
-  fc <- igraph::V(graph)$value
-  names(fc) <- igraph::V(graph)$name
-  if (!is.null(ntop) && ntop > 0) {
-    ii <- tapply(
-      seq_along(fc), igraph::V(graph)$layer,
-      function(i) utils::head(i[order(-abs(fc[i]))], ntop)
-    )
-    ii <- unlist(ii[names(ii) %in% layers])
-    fc <- fc[ii]
-    graph <- igraph::subgraph(graph, igraph::V(graph)[ii])
-  }
+  ## reduce the vertex set
+  if (!is.null(select)) graph <- keep_named_nodes(graph, select)
+  if (!is.null(filter)) graph <- keep_matching_nodes(graph, filter)
+  if (!is.null(ntop) && ntop > 0) graph <- keep_top_nodes(graph, ntop, layers)
+
+  ## zero the weight of every unwanted edge, then drop them in one go
+  if (normalize.edges) graph <- normalize_edge_weights(graph)
+  if (min.rho > 0) graph <- zero_weak_edges(graph, min.rho)
+  graph <- zero_edges_by_sign(graph, edge.sign, layers)
+  graph <- zero_edges_by_type(graph, edge.type)
+
+  return(drop_zero_edges(graph, prune))
+
+}
+
+## Layers the graph is to be reduced to: the caller's choice, otherwise the
+## layers attached to the graph object, otherwise those found on its
+## vertices. SOURCE and SINK are helper nodes and are never kept.
+resolve_layers <- function(graph, layers) {
+
+  if (is.null(layers)) layers <- graph$layers
+  if (is.null(layers)) layers <- unique(igraph::V(graph)$layer)
 
-  if (normalize.edges) {
-    for (e in unique(igraph::E(graph)$connection_type)) {
-      ii <- which(igraph::E(graph)$connection_type == e)
-      max.wt <- max(abs(igraph::E(graph)$weight[ii]), na.rm = TRUE) + 1e-3
-      igraph::E(graph)$weight[ii] <- igraph::E(graph)$weight[ii] / max.wt
-    }
+  return(setdiff(layers, c("SOURCE", "SINK")))
+
+}
+
+## Keep the nodes named in 'select', matching either the full node name or
+## the part before or after its layer prefix, so that a module can be
+## selected by bare name as well as by qualified name.
+keep_named_nodes <- function(graph, select) {
+
+  vv <- igraph::V(graph)$name
+  v1 <- (vv %in% select)
+  v2 <- (sub(".*:", "", vv) %in% select)
+  v3 <- (sub(":.*", "", vv) %in% select)
+
+  return(igraph::subgraph(graph, vids = which(v1 | v2 | v3)))
+
+}
+
+## Keep the nodes matching the per-layer regular expressions in 'filter', a
+## list of patterns named by layer. Each pattern only constrains its own
+## layer; nodes of the layers not named in the list are always kept.
+keep_matching_nodes <- function(graph, filter) {
+
+  for (k in names(filter)) {
+    vv <- igraph::V(graph)$name
+    filt <- filter[[k]]
+    vids <- igraph::V(graph)$layer != k | grepl(filt, vv, ignore.case = TRUE)
+    graph <- igraph::subgraph(graph, which(vids))
   }
 
-  if (min.rho > 0) {
-    ii <- which(abs(igraph::E(graph)$weight) < min.rho)
-    if (length(ii)) igraph::E(graph)$weight[ii] <- 0
+  return(graph)
+
+}
+
+## Keep the 'ntop' nodes of largest absolute value within each layer. Layers
+## outside 'layers' are dropped entirely, which matters for the SOURCE/SINK
+## nodes that survive an explicit layer selection.
+keep_top_nodes <- function(graph, ntop, layers) {
+
+  fc <- igraph::V(graph)$value
+  ii <- tapply(
+    seq_along(fc), igraph::V(graph)$layer,
+    function(i) utils::head(i[order(-abs(fc[i]))], ntop)
+  )
+  ii <- unlist(ii[names(ii) %in% layers])
+
+  return(igraph::subgraph(graph, igraph::V(graph)[ii]))
+
+}
+
+## Rescale the edge weights of each connection type onto a common scale, so
+## that weights of different connection types can be compared.
+normalize_edge_weights <- function(graph) {
+
+  for (e in unique(igraph::E(graph)$connection_type)) {
+    ii <- which(igraph::E(graph)$connection_type == e)
+    max.wt <- max(abs(igraph::E(graph)$weight[ii]), na.rm = TRUE) + 1e-3
+    igraph::E(graph)$weight[ii] <- igraph::E(graph)$weight[ii] / max.wt
   }
 
+  return(graph)
+
+}
+
+## Zero the edges weaker than 'min.rho'.
+zero_weak_edges <- function(graph, min.rho) {
+
+  ii <- which(abs(igraph::E(graph)$weight) < min.rho)
+  if (length(ii)) igraph::E(graph)$weight[ii] <- 0
+
+  return(graph)
+
+}
+
+## Zero the edges of unwanted sign. "consensus" keeps the edges whose sign
+## agrees with the expected direction of their source layer, which is
+## inverted for the miRNA layers as those anticorrelate with their targets.
+zero_edges_by_sign <- function(graph, edge.sign, layers) {
+
   ewt <- igraph::E(graph)$weight
+
   if (grepl("pos", edge.sign)) {
     igraph::E(graph)$weight[ewt < 0] <- 0
   } else if (grepl("neg", edge.sign)) {
@@ -111,8 +177,17 @@ prune_graph <- function(graph,
     igraph::E(graph)$weight <- ewt * (sign(ewt) == esign)
   }
 
-  ## delete intra or inter edges
+  return(graph)
+
+}
+
+## Zero the edges of unwanted type, inter-layer edges being those whose
+## connection type names two layers. "both2" keeps both types but drops the
+## negative intra-layer edges.
+zero_edges_by_type <- function(graph, edge.type) {
+
   ic <- grepl("->", igraph::E(graph)$connection_type)
+
   if (edge.type == "inter") {
     igraph::E(graph)$weight[!ic] <- 0
   } else if (edge.type == "intra") {
@@ -121,8 +196,16 @@ prune_graph <- function(graph,
     sel <- (!ic & igraph::E(graph)$weight < 0)
     igraph::E(graph)$weight[sel] <- 0
   }
-  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
 
+  return(graph)
+
+}
+
+## Drop the edges that the filters zeroed out and, with prune=TRUE, the
+## vertices that are left without any edge.
+drop_zero_edges <- function(graph, prune) {
+
+  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
   if (prune) graph <- igraph::subgraph_from_edges(graph, igraph::E(graph))
 
   return(graph)
diff --git a/R/lasagna-package.R b/R/lasagna-package.R
index 443dd6b..643d87e 100644
--- a/R/lasagna-package.R
+++ b/R/lasagna-package.R
@@ -2,6 +2,7 @@
 "_PACKAGE"
 
 #' @import igraph
+#' @importFrom magrittr %>%
 #' @importFrom matrixStats rowSds
 #' @importFrom stats cor median model.matrix sd
 #' @importFrom utils head tail type.convert
diff --git a/R/model.R b/R/model.R
index 0007cdd..1d60c5b 100644
--- a/R/model.R
+++ b/R/model.R
@@ -6,6 +6,11 @@
 #' Edges weighted by correlation, optionally conditioned on phenotype.
 #' @param data A list with \code{X} (named list of data matrices),
 #'   \code{samples} (data frame), and optionally \code{contrasts}.
+#' @param X Named list of data matrices, one per layer. Alternative to
+#'   \code{data} when \code{data} is not supplied.
+#' @param meta Sample phenotype data frame (or contrasts matrix, per
+#'   \code{meta.type}). Alternative to \code{data} when \code{data} is
+#'   not supplied.
 #' @param meta.type Phenotype type: \code{"pheno"}, \code{"expanded"},
 #'   or \code{"contrasts"}.
 #' @param ntop Number of top-SD features per layer. Set 0 or NULL
@@ -53,7 +58,7 @@ create_model <- function(data,
                          fully_connect = FALSE,
                          add.revpheno = TRUE,
                          condition.edges = TRUE) {
-  
+
   if(!is.null(data)) {
     X <- data$X
     meta <- data$samples
@@ -68,9 +73,13 @@ create_model <- function(data,
   if(is.null(data) && (is.null(X) || is.null(meta)) ) {
     stop("must supply data or {X, meta}.")
   }
-  
+
   if (meta.type %in% c("pheno","samples")) {
     Y <- expandPhenoMatrix(meta, drop.ref = FALSE)
+    if (is.null(Y)) {
+      stop("[create_model] no phenotype column with resolvable/varying groups found in samples/meta; ",
+           "check that it has a factor or character column with at least 2 distinct values")
+    }
   } else if (meta.type %in% c("expanded","traits")) {
     Y <- 1 * meta
   } else if (meta.type == "contrasts") {
@@ -93,12 +102,35 @@ create_model <- function(data,
   xx <- X
   xx <- lapply(xx, as.matrix)
   if (!is.null(ntop) && ntop > 0) {
-    xx <- lapply(xx, function(x) head(x[order(-apply(x, 1, stats::sd)), , drop = FALSE], ntop))
     xx <- mofa.topSD(xx, ntop)
   }
 
   ## merge data (handles non-matching samples)
   xx <- mofa.merge_data2(xx, merge.rows = "prefix", merge.cols = "union")
+
+  if (ncol(xx)<2) {
+    message("[create_model] < 2 samples detected. Cannot compute correlation. Returning NULL. Exiting.")
+    return(NULL)
+  }
+
+  ## diagnose sample overlap across layers: warn (don't fail) when most
+  ## samples lack data in at least one layer, since this silently degrades
+  ## downstream correlation/model fit
+  row.layer <- sub(":.*", "", rownames(xx))
+  layer.names <- unique(row.layer)
+  layer.coverage <- sapply(layer.names, function(lyr) {
+    rows <- which(row.layer == lyr)
+    colSums(!is.na(xx[rows, , drop = FALSE])) > 0
+  })
+  complete.frac <- mean(rowSums(layer.coverage) == length(layer.names))
+  na.frac <- mean(is.na(xx))
+  if (complete.frac < 0.5) {
+    warning(sprintf(
+      "create_model: merged data is %.0f%% NA after combining %d layers with only partial sample overlap; only %.0f%% of samples have data in every layer; results may be unreliable",
+      100 * na.frac, length(layer.names), 100 * complete.frac
+    ), call. = FALSE)
+  }
+
   kk <- intersect(colnames(xx), rownames(Y))
   xx <- xx[, kk]
   Y <- Y[kk, ]
@@ -123,7 +155,8 @@ create_model <- function(data,
   if (condition.edges) {
     message("conditioning edges...")
     rho <- stats::cor(t(xx), Y, use = "pairwise.complete.obs")
-    maxrho <- apply(abs(rho), 1, max, na.rm = TRUE)
+    maxrho <- matrixStats::rowMaxs(abs(rho), na.rm = TRUE)
+    names(maxrho) <- rownames(rho)
     ii <- grep("SINK|SOURCE", names(maxrho))
     if (length(ii)) maxrho[ii] <- 1
     rho.wt <- outer(maxrho, maxrho)
@@ -139,11 +172,13 @@ create_model <- function(data,
   if (!fully_connect) {
     layer_mask <- matrix(0, nrow(R), ncol(R))
     dimnames(layer_mask) <- dimnames(R)
-    for (i in seq_len(length(layers) - 1)) {
-      ii <- which(dt == layers[i])
-      jj <- which(dt == layers[i + 1])
-      layer_mask[ii, jj] <- 1
-      layer_mask[jj, ii] <- 1
+    if (length(layers) >= 2) {
+      for (i in seq_len(length(layers) - 1)) {
+        ii <- which(dt == layers[i])
+        jj <- which(dt == layers[i + 1])
+        layer_mask[ii, jj] <- 1
+        layer_mask[jj, ii] <- 1
+      }
     }
     if (intra) {
       for (i in seq_along(layers)) {
@@ -159,22 +194,28 @@ create_model <- function(data,
     message("reducing edges to maximum ", nc, " connections")
     xtypes <- setdiff(layers, c("PHENO", "SOURCE", "SINK"))
     reduce_mask <- matrix(1, nrow(R), ncol(R))
-    for (i in seq_len(length(xtypes) - 1)) {
-      ii <- which(dt == xtypes[i])
-      jj <- which(dt == xtypes[i + 1])
-      R1 <- R[ii, jj, drop = FALSE]
-      rii <- apply(abs(R1), 1, function(r) utils::tail(sort(r), nc)[1])
-      rjj <- apply(abs(R1), 2, function(r) utils::tail(sort(r), nc)[1])
-      rr <- abs(R1) >= rii | t(t(abs(R1)) >= rjj)
-      reduce_mask[ii, jj] <- rr
-      reduce_mask[jj, ii] <- t(rr)
+    if (length(xtypes) >= 2) {
+      for (i in seq_len(length(xtypes) - 1)) {
+        ii <- which(dt == xtypes[i])
+        jj <- which(dt == xtypes[i + 1])
+        R1 <- R[ii, jj, drop = FALSE]
+        k_row <- pmax(ncol(R1) - nc + 1, 1)
+        k_col <- pmax(nrow(R1) - nc + 1, 1)
+        rii <- matrixStats::rowOrderStats(abs(R1), which = k_row)
+        rjj <- matrixStats::colOrderStats(abs(R1), which = k_col)
+        rr <- abs(R1) >= rii | t(t(abs(R1)) >= rjj)
+        reduce_mask[ii, jj] <- rr
+        reduce_mask[jj, ii] <- t(rr)
+      }
     }
     if (intra) {
       for (i in seq_along(xtypes)) {
         ii <- which(dt == xtypes[i])
         R1 <- R[ii, ii, drop = FALSE]
-        rii <- apply(abs(R1), 1, function(r) utils::tail(sort(r), nc)[1])
-        rjj <- apply(abs(R1), 2, function(r) utils::tail(sort(r), nc)[1])
+        k_row <- pmax(ncol(R1) - nc + 1, 1)
+        k_col <- pmax(nrow(R1) - nc + 1, 1)
+        rii <- matrixStats::rowOrderStats(abs(R1), which = k_row)
+        rjj <- matrixStats::colOrderStats(abs(R1), which = k_col)
         rr <- abs(R1) >= rii | t(t(abs(R1)) >= rjj)
         reduce_mask[ii, ii] <- rr
       }
@@ -194,16 +235,26 @@ create_model <- function(data,
   ## add edge connection type as attribute
   igraph::V(gr)$layer <- sub(":.*", "", igraph::V(gr)$name)
   ee <- igraph::as_edgelist(gr)
-  etype <- apply(ee, 2, function(e) sub(":.*", "", e))
-  etype.idx <- apply(etype, 2, match, gr$layers)
-  rev.etype <- etype.idx[, 2] < etype.idx[, 1]
-  etype1 <- ifelse(rev.etype, etype.idx[, 2], etype.idx[, 1])
-  etype2 <- ifelse(rev.etype, etype.idx[, 1], etype.idx[, 2])
-  etype1 <- gr$layers[etype1]
-  etype2 <- gr$layers[etype2]
-  igraph::E(gr)$connection_type <- paste0(etype1, "->", etype2)
-  ii <- which(etype1 == etype2)
-  if (length(ii)) igraph::E(gr)$connection_type[ii] <- etype1[ii]
+  if (nrow(ee) < 2) {
+    ## apply() collapses to a plain vector (losing the edge x endpoint
+    ## matrix shape) when there are 0 or 1 edges; handle those directly
+    igraph::E(gr)$connection_type <- vapply(seq_len(nrow(ee)), function(i) {
+      ty <- sort(match(sub(":.*", "", ee[i, ]), gr$layers))
+      lyr <- gr$layers[ty]
+      if (lyr[1] == lyr[2]) lyr[1] else paste0(lyr[1], "->", lyr[2])
+    }, character(1))
+  } else {
+    etype <- apply(ee, 2, function(e) sub(":.*", "", e))
+    etype.idx <- apply(etype, 2, match, gr$layers)
+    rev.etype <- etype.idx[, 2] < etype.idx[, 1]
+    etype1 <- ifelse(rev.etype, etype.idx[, 2], etype.idx[, 1])
+    etype2 <- ifelse(rev.etype, etype.idx[, 1], etype.idx[, 2])
+    etype1 <- gr$layers[etype1]
+    etype2 <- gr$layers[etype2]
+    igraph::E(gr)$connection_type <- paste0(etype1, "->", etype2)
+    ii <- which(etype1 == etype2)
+    if (length(ii)) igraph::E(gr)$connection_type[ii] <- etype1[ii]
+  }
 
   return(list(graph = gr, X = xx, Y = Y, layers = layers))
 
diff --git a/R/plot.R b/R/plot.R
index 73902c8..1dc66bf 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -17,6 +17,7 @@
 #' @param justgraph Only plot graph (no labels/titles).
 #' @param edge.cex Edge width multiplier.
 #' @param edge.alpha Edge transparency.
+#' @param edge.colors Two colors for negative and positive edge weights.
 #' @param xdist Distance between layers.
 #' @param normalize.edges Normalize edges per connection type.
 #' @param yheight Height of the y-axis layout.
@@ -75,6 +76,11 @@ plot_multipartite <- function(graph,
                               color.var = "value",
                               layout = c("parallel", "hive")[1]) {
 
+  if (is.null(graph)) {
+    message("[lasagna::plot_multipartite]: graph is NULL. Please check run of create_model() and solve(). Exiting.")
+    return(NULL)
+  }
+
   vattr <- igraph::vertex_attr_names(graph)
   edgeattr <- igraph::edge_attr_names(graph)
   if (!"rho" %in% edgeattr) warning("no rho in edge attributes")
@@ -96,8 +102,9 @@ plot_multipartite <- function(graph,
     prune = prune
   )
 
-  if (length(igraph::V(graph)) == 0) warning("graph has no nodes")
-  if (length(igraph::E(graph)) == 0) warning("graph has no edges")
+  if (length(igraph::V(graph)) == 0 || length(igraph::E(graph)) == 0) {
+    stop("plot_multipartite: no nodes/edges remain after filtering (layers=/ntop=/min.rho=) -- nothing to plot")
+  }
   
   layers <- graph$layers
   layers <- setdiff(layers, c("SOURCE", "SINK"))
@@ -244,8 +251,12 @@ plot_multipartite <- function(graph,
 #' @param egamma Gamma exponent applied to edge widths.
 #' @param color.var Vertex colouring: \code{"value"}, \code{"type"}/\code{"layer"},
 #'   or \code{"color"} (use the vertex \code{color} attribute).
+#' @param edge.colors Two colors for negative and positive edge weights.
 #' @param labcex Label size multiplier.
-#' @param layout Optional layout matrix (rows named by vertex).
+#' @param layout Optional layout matrix (rows named by vertex). If
+#'   \code{NULL}, a \code{layout} attribute pre-set on \code{graph}
+#'   (\code{graph$layout <- ...}) is used before falling back to the
+#'   default visNetwork layout.
 #' @param physics Enable physics simulation.
 #' @return A visNetwork widget.
 #' @examples
@@ -275,6 +286,74 @@ plot_visgraph <- function(graph,
                           layout = NULL,
                           physics = TRUE) {
 
+  if (is.null(graph)) {
+    message("[lasagna::plot_visgraph]: graph is NULL. Please check run of create_model() and solve(). Nothing to plot. Exiting.")
+    return(NULL)
+  }
+
+  sub <- visgraph_subgraph(graph, layers, ntop, min_rho, mst)
+  sub <- visgraph_style(sub, color.var, edge.colors, labcex, ecex, egamma)
+
+  data <- visNetwork::toVisNetworkData(sub, idToLabel=FALSE)
+
+  vis <- visNetwork::visNetwork(
+    nodes = data$nodes,
+    edges = data$edges,
+    height = "800px", width = "100%"
+  ) %>%
+    visNetwork::visNodes(
+      scaling = list(min = 5 * vcex, max = 15*vcex),
+      shadow = list(enable =TRUE),      
+      color = list(border = "black"),
+      font = list(align = "left")
+    ) %>%
+    visNetwork::visEdges(
+      color = list(opacity = 0.2),
+      scaling = list(min = 3, max = 30)
+    ) %>%
+    visNetwork::visInteraction(
+      hover = TRUE
+    ) %>%
+    visNetwork::visOptions(
+      highlightNearest = TRUE
+    ) %>%
+    visNetwork::visPhysics(
+      enable = physics,
+      barnesHut = list(
+        springLength = 50
+      )
+    )
+
+  M <- visgraph_layout_matrix(layout, sub, data$nodes$id)
+
+  if (!is.null(M)) {
+    vis <- vis %>% visNetwork::visIgraphLayout(layout = "layout.norm", layoutMatrix = M)
+  }
+
+  return(vis)
+
+}
+
+## Coordinate matrix to lay the nodes out with, in node order: the caller's
+## layout, otherwise a layout pre-attached to the graph object, with y
+## flipped to match the visNetwork axis. NULL when neither is given.
+visgraph_layout_matrix <- function(layout, sub, ids) {
+
+  if (is.null(layout)) layout <- sub$layout
+  if (is.null(layout)) return(NULL)
+
+  M <- layout[ids,]
+  M[,2] <- -M[,2]
+
+  return(M)
+
+}
+
+## Reduce the graph to what plot_visgraph will actually draw: the requested
+## layers, optionally its minimum spanning tree, the ntop strongest
+## vertices, and the edges above min_rho.
+visgraph_subgraph <- function(graph, layers, ntop, min_rho, mst) {
+
   if (is.null(layers)) layers <- graph$layers
 
   sub <- igraph::subgraph(graph, igraph::V(graph)$layer %in% layers)
@@ -293,6 +372,15 @@ plot_visgraph <- function(graph,
     sub <- igraph::subgraph_from_edges(sub, which(abs(igraph::E(sub)$weight) > min_rho))
   }
 
+  return(sub)
+
+}
+
+## Set the visNetwork drawing attributes on the subgraph: node color by the
+## chosen variable, up/down triangle shapes, node sizes from the (squared)
+## node values, and edge width/color from the edge weights.
+visgraph_style <- function(sub, color.var, edge.colors, labcex, ecex, egamma) {
+
   vtype <- sub(":.*", "", igraph::V(sub)$name)
   ntypes <- length(unique(vtype))
   vcol <- "grey"
@@ -311,7 +399,7 @@ plot_visgraph <- function(graph,
   vcol <- sub("^[A-Z]+","",vcol) ## strip away WGCNA prefix
   igraph::V(sub)$color <- vcol
   igraph::V(sub)$color.border <- "black"
- 
+
   vtype <- c("down", "up")[1 + 1 * (igraph::V(sub)$value > 0)]
   igraph::V(sub)$shape <- c("triangleDown", "triangle")[as.factor(vtype)]
   igraph::V(sub)$value <- abs(igraph::V(sub)$value)^2
@@ -326,47 +414,8 @@ plot_visgraph <- function(graph,
 
   igraph::E(sub)$width <- 5 * ecex * abs(igraph::E(sub)$weight)^egamma
   igraph::E(sub)$color <- edge.colors[1 + 1 * (igraph::E(sub)$weight > 0)]
-  
-  data <- visNetwork::toVisNetworkData(sub, idToLabel=FALSE)
 
-  vis <- visNetwork::visNetwork(
-    nodes = data$nodes,
-    edges = data$edges,
-    height = "800px", width = "100%"
-  ) %>%
-    visNetwork::visNodes(
-      scaling = list(min = 5 * vcex, max = 15*vcex),
-      shadow = list(enable =TRUE),      
-      color = list(border = "black"),
-      font = list(align = "left")
-    ) %>%
-    visNetwork::visEdges(
-      color = list(opacity = 0.2),
-      scaling = list(min = 3, max = 30)
-    ) %>%
-    visNetwork::visInteraction(
-      hover = TRUE
-    ) %>%
-    visNetwork::visOptions(
-      highlightNearest = TRUE
-    ) %>%
-    visNetwork::visPhysics(
-      enable = physics,
-      barnesHut = list(
-        springLength = 50
-      )
-    )
-
-  if (is.null(layout) && !is.null(sub$layout)) layout <- sub$layout
-  
-  if (!is.null(layout)) {
-    vv <- data$nodes$id
-    M <- layout[vv,]
-    M[,2] <- -M[,2]
-    vis <- vis %>% visNetwork::visIgraphLayout(layout = "layout.norm", layoutMatrix = M)
-  }
-  
-  return(vis)
+  return(sub)
 
 }
 
@@ -377,7 +426,12 @@ plot_visgraph <- function(graph,
 #' @param graph An igraph object (output of \code{solve}).
 #' @param layout Either a matrix/data frame with columns \code{x}, \code{y},
 #'   \code{z} (one row per vertex, row names matching vertex names), or a named
-#'   list of 2-column position matrices per layer.
+#'   list of 2-column position matrices per layer. If \code{NULL}, a
+#'   \code{layout} attribute pre-set on \code{graph} (\code{graph$layout <-
+#'   ...}) is used before falling back to an automatically computed layout.
+#' @param X Numeric matrix of features (rows, named to match graph vertex
+#'   names) by samples (columns), used to compute the layout when
+#'   \code{layout} is \code{NULL} or a method name.
 #' @param draw_edges Logical; draw inter-layer edges.
 #' @param num_edges Maximum number of edges per layer pair.
 #' @param min_rho Minimum absolute weight for edges.
@@ -385,8 +439,12 @@ plot_visgraph <- function(graph,
 #'   \code{"both"}.
 #' @param cex Point size multiplier.
 #' @param cex.gamma Gamma exponent applied to point sizes.
+#' @param edge.cex Edge width multiplier.
 #' @param color.by Vertex colouring: \code{"value"} or \code{"color"}.
+#' @param edge_colors Two colors for negative and positive edge weights.
 #' @param znames Named character vector mapping layer codes to display names.
+#' @param ax Integer (mod 4) selecting which plot corner the layer title
+#'   text is placed at.
 #' @return A plotly object.
 #' @examples
 #' set.seed(1)
@@ -420,37 +478,78 @@ plot_3d <- function(graph,
                     color.by = "value",
                     edge_colors = c("blue", "magenta"),
                     znames = NULL, ax=0) {
-  require(plotly)
-  
+
+  if (is.null(graph)) {
+    message("[lasagna::plot_3d]: graph is NULL. Please check run of create_model() and solve(). Nothing to plot. Exiting.")
+    return(NULL)
+  }
+
   edges <- NULL
   if (draw_edges) {
     edges <- data.frame(igraph::as_edgelist(graph),
       weight = igraph::E(graph)$weight)
   }
 
+  layout <- resolve_layout_3d(graph, layout, X)
+
+  ## remove SINK/SOURCE
+  vv <- intersect(c("SINK","SOURCE"), igraph::V(graph)$name)
+  if(length(vv)) graph <- igraph::delete_vertices(graph, vv)
+
+  df <- build_node_frame_3d(graph, layout, cex.gamma, color.by)
+  edges <- select_edges_3d(edges, min_rho, sign_rho, num_edges)
+  if (is.null(znames)) znames <- default_layer_names()
+
+  edge.colors <- c()
+
+  plt <- plotlyLasagna(df, znames = znames, edges = edges,
+    edge_colors = edge_colors,
+    cex = cex, edge.cex = edge.cex, ax = ax)
+
+  return(plt)
+}
+
+## Resolve the 'layout' argument of plot_3d into a coordinate table: a
+## user-supplied data frame, a layout pre-attached to the graph object, or
+## a layout computed from X by the requested reduction method. These are
+## checks on plot_3d's own arguments, so they are reported against the
+## caller rather than against this helper.
+resolve_layout_3d <- function(graph, layout, X) {
+
+  caller <- sys.call(-1)
+  abort <- function(msg) stop(simpleError(msg, call = caller))
+
   vx <- igraph::V(graph)$name
+
   if(!is.null(layout) && is.data.frame(layout)) {
     message("using user layout")
-  } else if(!is.null(layout) && is.character(layout)) {
-    if(is.null(X)) stop("must provide X or layout")
-    if(!all(vx %in% rownames(X))) stop("incomplete layout")    
-    if(!layout %in% c("svd","tsne","umap")) stop("invalid layout")
-    layout <- layout_multipartite_3d(graph, X, clust=layout)
-  } else if(is.null(layout) && !is.null(graph$layout))  {
-    message("using layout in graph object")
-    layout <- graph$layout
-  } else {
-    if(is.null(X)) stop("must provide X or layout")
-    if(!all(vx %in% rownames(X))) stop("incomplete layout")
-    layout <- layout_multipartite_3d(graph, X, clust="umap")
+    return(layout)
   }
 
-  ## remove SINK/SOURCE
-  if(1) {
-    vv <- intersect(c("SINK","SOURCE"), igraph::V(graph)$name)
-    if(length(vv)) graph <- igraph::delete_vertices(graph, vv)
+  if(!is.null(layout) && is.character(layout)) {
+    if(is.null(X)) abort("must provide X or layout")
+    if(!all(vx %in% rownames(X))) abort("incomplete layout")
+    if(!layout %in% c("svd","tsne","umap")) abort("invalid layout")
+    return(layout_multipartite_3d(graph, X, clust=layout))
+  }
+
+  if(is.null(layout) && !is.null(graph$layout)) {
+    ## allow a layout pre-attached to the graph object to override the default
+    message("using layout in graph object")
+    return(graph$layout)
   }
 
+  if(is.null(X)) abort("must provide X or layout")
+  if(!all(vx %in% rownames(X))) abort("incomplete layout")
+
+  return(layout_multipartite_3d(graph, X, clust="umap"))
+
+}
+
+## Assemble the per-vertex plotting frame for plot_3d: layer-wise rescaled
+## x/y positions plus value, point size, color and hover text.
+build_node_frame_3d <- function(graph, layout, cex.gamma, color.by) {
+
   ## feature maps across datatypes
   vx <- igraph::V(graph)$name
   df <- data.frame(layout[vx,])
@@ -469,7 +568,7 @@ plot_3d <- function(graph,
   df$value <- as.numeric(vars)
   df$size <- abs(as.numeric(vars))^cex.gamma
   df$color <- as.numeric(vars)
-  
+
   if (!is.null(igraph::V(graph)$color) && color.by == "color") {
     df$color <- igraph::V(graph)$color
     df$color <- sub("^[A-Z]+","",df$color) ## remove prefix
@@ -484,8 +583,17 @@ plot_3d <- function(graph,
   #df$z <- factor(df$z, levels = levels)
   df$text <- paste(rownames(df), "<br>value:", round(df$value, digits = 3))
 
-  ## filter edges
-  if (!is.null(edges) && min_rho >= 0) {
+  return(df)
+
+}
+
+## Keep only the edges that plot_3d will draw: filter on sign and minimum
+## weight, then cap the number of edges per layer pair.
+select_edges_3d <- function(edges, min_rho, sign_rho, num_edges) {
+
+  if (is.null(edges)) return(NULL)
+
+  if (min_rho >= 0) {
     if (sign_rho == "pos") {
       edges <- edges[ which(edges[, 3] > min_rho), ]
     } else if (sign_rho == "neg") {
@@ -495,7 +603,7 @@ plot_3d <- function(graph,
     }
   }
 
-  if (!is.null(edges) && num_edges > 0) {
+  if (num_edges > 0) {
     e1 <- sub(":.*","",edges[,1])
     e2 <- sub(":.*","",edges[,2])
     ee <- paste0(e1,'-',e2)
@@ -508,43 +616,42 @@ plot_3d <- function(graph,
     edges <- edges[edges[, 3] != 0, ]
   }
 
-  ## layer name mapping
-  if (is.null(znames)) {
-    znames <- c(
-      "PHENO" = "Phenotype",
-      "ph" = "Phenotype",
-      "gset" = "Pathway",
-      "mx" = "Metabolomics",
-      "lx" = "Lipidomics",      
-      "gx" = "Transcriptomics",
-      "tx" = "Transcriptomics",
-      "mir" = "microRNA",
-      "px" = "Proteomics",
-      "hx" = "Histone",
-      "hptm" = "hPTM",
-      "dr" = "Drug response",
-      "me" = "Methylation",
-      "mt" = "Mutation",
-      "mut" = "Mutation",      
-      "mu" = "Mutation"
-    )
-  }
+  return(edges)
 
-  edge.colors <- c()
-  
-  plt <- plotlyLasagna(df, znames = znames, edges = edges,
-    edge_colors = edge_colors,
-    cex = cex, edge.cex = edge.cex, ax = ax)
-  
-  return(plt)
+}
+
+## Default mapping of layer codes to display names.
+default_layer_names <- function() {
+  c(
+    "PHENO" = "Phenotype",
+    "ph" = "Phenotype",
+    "gset" = "Pathway",
+    "mx" = "Metabolomics",
+    "lx" = "Lipidomics",
+    "gx" = "Transcriptomics",
+    "tx" = "Transcriptomics",
+    "mir" = "microRNA",
+    "px" = "Proteomics",
+    "hx" = "Histone",
+    "hptm" = "hPTM",
+    "dr" = "Drug response",
+    "me" = "Methylation",
+    "mt" = "Mutation",
+    "mut" = "Mutation",
+    "mu" = "Mutation"
+  )
 }
 
 #' Internal plotly builder for LASAGNA 3D plot
 #' Builds the actual plotly figure from a prepared data frame.
 #' @param df Data frame with columns: feature, x, y, z, color, text.
 #' @param znames Named character vector for layer display names.
+#' @param ax Integer (mod 4) selecting which plot corner the layer title
+#'   text is placed at.
 #' @param cex Point size multiplier.
+#' @param edge.cex Edge line width multiplier.
 #' @param edges Optional data frame of edges.
+#' @param edge_colors Two colors for negative and positive edge weights.
 #' @return A plotly object.
 #' @examples
 #' set.seed(1)
@@ -560,9 +667,12 @@ plotlyLasagna <- function(df,
                           znames = NULL, ax = 1,
                           cex = 1, edge.cex = 1,
                           edges = NULL,
-                          edge_colors = c("blue", "magenta")
-                          ) {
-  require(plotly)
+                          edge_colors = c("blue", "magenta")) {
+
+  if (is.null(df) || nrow(df) == 0) {
+    message("[lasagna::plotlyLasagna]: df is NULL or empty. Nothing to plot. Exiting.")
+    return(NULL)
+  }
 
   zz <- sort(unique(df$z))
   min.x <- min(df$x, na.rm = TRUE)
@@ -623,18 +733,10 @@ plotlyLasagna <- function(df,
 
     ## add segments
     if (k < length(zz) && !is.null(edges)) {
-      sel1 <- which(edgetype1 == zz[k] & edgetype2 == zz[k + 1])
-      sel2 <- which(edgetype2 == zz[k] & edgetype1 == zz[k + 1])
-      df2 <- df[which(df$z == zz[k + 1]), c("x", "y", "z")]
-      sel <- unique(c(sel1, sel2))
-
-      if (length(sel)) {
-        ee <- edges[sel, ]
-        idx <- as.vector(t(as.matrix(ee[, 1:2])))
-        dfe <- rbind(df1[, c("x", "y", "z")], df2[, c("x", "y", "z")])[idx, ]
-        dfe$pair_id <- as.vector(mapply(rep, 1:nrow(ee), 2))
-        cc <- edge_colors[1 + (ee[, 3] > 0)]
-        dfe$col <- as.vector(mapply(rep, cc, 2))
+      dfe <- lasagna_edge_segments(df, edges, edgetype1, edgetype2,
+        zz[k], zz[k + 1], edge_colors)
+
+      if (!is.null(dfe)) {
         fig <- fig %>%
           plotly::add_trace(
             x = dfe$x,
@@ -651,25 +753,16 @@ plotlyLasagna <- function(df,
     }
 
     ## add layer title
-    ztext <- z
-    if (!is.null(znames) && (is.integer(z) || z %in% names(znames))) {
-      if (is.factor(z)) z <- as.character(z)
-      ztext <- znames[z]
-    }
+    title <- lasagna_layer_title(z, znames)
+    corner <- lasagna_axis_corner(ax, min.x, max.x, min.y, max.y)
 
-    ax <- ax %% 4
-    if(ax==0) {zx=min.x; zy=max.y}
-    if(ax==1) {zx=min.x; zy=min.y}
-    if(ax==2) {zx=max.x; zy=max.y}
-    if(ax==3) {zx=max.x; zy=min.y}
-    
     fig <- fig %>%
       plotly::add_text(
-        x = zx,
-        y = zy,
-        z = z,
+        x = corner$x,
+        y = corner$y,
+        z = title$z,
         mode = "text",
-        text = ztext,
+        text = title$text,
         textfont = list(size = 24),
         showlegend = FALSE,
         inherit = FALSE
@@ -690,6 +783,61 @@ plotlyLasagna <- function(df,
   fig
 }
 
+## Build the line segments joining two adjacent layers: the edges running
+## between them, laid out as one x/y/z row per endpoint, grouped by pair
+## and colored by the sign of the edge weight. Returns NULL if the two
+## layers share no edges.
+lasagna_edge_segments <- function(df, edges, edgetype1, edgetype2, z1, z2,
+                                  edge_colors) {
+
+  sel1 <- which(edgetype1 == z1 & edgetype2 == z2)
+  sel2 <- which(edgetype2 == z1 & edgetype1 == z2)
+  sel <- unique(c(sel1, sel2))
+  if (!length(sel)) return(NULL)
+
+  df1 <- df[which(df$z == z1), c("x", "y", "z")]
+  df2 <- df[which(df$z == z2), c("x", "y", "z")]
+
+  ee <- edges[sel, ]
+  idx <- as.vector(t(as.matrix(ee[, 1:2])))
+  dfe <- rbind(df1, df2)[idx, ]
+  dfe$pair_id <- rep(seq_len(nrow(ee)), each = 2)
+  cc <- edge_colors[1 + (ee[, 3] > 0)]
+  dfe$col <- rep(cc, each = 2)
+
+  return(dfe)
+
+}
+
+## Display name of a layer, and the z-coordinate to place it at. Looking up
+## a layer in znames also converts a factor level to its character label,
+## which the caller uses as the z-coordinate of the title.
+lasagna_layer_title <- function(z, znames) {
+
+  ztext <- z
+  if (!is.null(znames) && (is.integer(z) || z %in% names(znames))) {
+    if (is.factor(z)) z <- as.character(z)
+    ztext <- znames[z]
+  }
+
+  return(list(z = z, text = ztext))
+
+}
+
+## Plot corner (in x/y) that layer titles are anchored to, selected by
+## 'ax' modulo the four corners of the layer plane.
+lasagna_axis_corner <- function(ax, min.x, max.x, min.y, max.y) {
+
+  ax <- ax %% 4
+  if(ax==0) {zx=min.x; zy=max.y}
+  if(ax==1) {zx=min.x; zy=min.y}
+  if(ax==2) {zx=max.x; zy=max.y}
+  if(ax==3) {zx=max.x; zy=min.y}
+
+  return(list(x = zx, y = zy))
+
+}
+
 layout_multipartite <- function(graph, xpos=NULL, xdist=1) {
 
   layers <- graph$layers
@@ -739,8 +887,27 @@ layout_hiveplot <- function(graph) {
 
 }
 
+#' Compute 3D layout coordinates per layer via dimensionality reduction
+#' Splits the input data matrix by layer and applies SVD, t-SNE, or UMAP
+#' to compute 2D (x, y) coordinates within each layer, combined with a
+#' layer factor as the z-coordinate, for use as a layout by \code{plot_3d}.
+#' @param graph An igraph object (output of \code{solve}).
+#' @param X Numeric matrix of features (rows, named to match graph vertex
+#'   names) by samples (columns), spanning all layers.
+#' @param clust Dimensionality reduction method used within each layer:
+#'   \code{"svd"}, \code{"tsne"}, or \code{"umap"}.
+#' @return A data frame with columns \code{x}, \code{y} (per-layer 2D
+#'   layout coordinates) and \code{z} (a factor giving the layer of each
+#'   vertex), one row per graph vertex, row names matching vertex names.
 #' @export
 layout_multipartite_3d <- function(graph, X, clust=c("svd","tsne","umap")) {
+
+  if (is.null(graph) || igraph::vcount(graph) == 0) {
+    message("[lasagna::layout_multipartite_3d]: graph is NULL or empty. Nothing to plot. Exiting.")
+    return(NULL)
+  }
+
+  clust <- match.arg(clust)
   layers <- graph$layers
   layers <- setdiff(layers, c("SOURCE","SINK"))
 
@@ -751,11 +918,15 @@ layout_multipartite_3d <- function(graph, X, clust=c("svd","tsne","umap")) {
   for(k in names(ff)) {
     nn <- nrow(ff[[k]])
     if(nn > 5 && clust == 'tsne') {
-      px <- max(min(30, nn/3),2)
+      if (!requireNamespace("Rtsne", quietly = TRUE)) stop("package Rtsne required for clust='tsne'")
+      px <- max(min(30, (nn-1)/3),2)
       xy[[k]] <- Rtsne::Rtsne(ff[[k]], perplexity=px)$Y
     } else if(nn > 5 && clust == 'umap') {
+      if (!requireNamespace("uwot", quietly = TRUE)) stop("package uwot required for clust='umap'")
       nb <- max(min(15, nn/3),2)
       xy[[k]] <- uwot::umap(ff[[k]], n_neighbors=nb)
+    } else if(nn == 1) {
+      xy[[k]] <- matrix(0, nrow=1, ncol=2)
     } else {
       xy[[k]] <- svd(ff[[k]])$u[,1:2]
     }
diff --git a/R/solve.R b/R/solve.R
index 330753c..7261189 100644
--- a/R/solve.R
+++ b/R/solve.R
@@ -39,6 +39,11 @@ solve <- function(obj,
                   sp.weight = FALSE,
                   graph = NULL) {
 
+  if (is.null(obj)) {
+    message("[lasagna::solve]: obj is NULL. Please check run of create_model(). Exiting.")
+    return(NULL)
+  }
+
   if (!pheno %in% colnames(obj$Y)) stop("pheno not in Y")
 
   if (!"rho" %in% names(igraph::edge_attr(obj$graph))) {
@@ -47,10 +52,37 @@ solve <- function(obj,
 
   if (is.null(graph)) graph <- obj$graph
 
-  X <- obj$X
-  y <- obj$Y[, pheno]
+  ## score the nodes against the phenotype
+  dat <- mask_neutral_samples(obj$X, obj$Y[, pheno])
+  rho <- node_correlation(dat$X, dat$y)
+  fc <- node_foldchange(dat$X, dat$y, rho, pheno)
+
+  ## weight the edges from the node scores
+  graph <- set_node_values(graph, rho, fc, value.type)
+  graph <- set_edge_weights(graph, fc.weights)
+  if (sp.weight) graph <- add_sp_weight(graph, obj$layers)
+
+  ## zero the weight of every unwanted edge, then drop them in one go
+  if (min_rho > 0) graph <- zero_weak_edges(graph, min_rho)
+  if (max_edges > 0) graph <- zero_excess_edges(graph, max_edges)
+  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
+
+  ## prune vertices
+  if (prune) {
+    ewt <- igraph::E(graph)$weight
+    graph <- igraph::subgraph_from_edges(graph, which(abs(ewt) > 0))
+  }
+
+  return(graph)
+
+}
+
+## A phenotype coded -1/0/1 uses the zeros for the samples that belong to
+## neither side of the contrast. Mask those out of both the phenotype and its
+## PHENO rows in the data, so that they weigh on neither correlation nor fold
+## change.
+mask_neutral_samples <- function(X, y) {
 
-  ## check if phenotype was coded -1/0/1
   ii <- grep("PHENO", rownames(X))
   has.min1 <- (min(X[ii, ], na.rm = TRUE) < 0)
   if (length(ii) && has.min1) {
@@ -58,29 +90,61 @@ solve <- function(obj,
     y[y == 0] <- NA
   }
 
+  return(list(X = X, y = y))
+
+}
+
+## Correlation of each node with the phenotype. Nodes that are constant or
+## never observed together with the phenotype score zero rather than NA.
+node_correlation <- function(X, y) {
+
   rho <- stats::cor(t(X), y, use = "pairwise")[, 1]
+  rho[is.na(rho)] <- 0
+
+  return(rho)
+
+}
+
+## Fold change of each node between the two sides of the phenotype. For the
+## PHENO nodes a fold change is meaningless, so they take their correlation
+## instead.
+node_foldchange <- function(X, y, rho, pheno) {
+
   i0 <- which(y <= 0)
   i1 <- which(y > 0)
+  if (length(i0) == 0 || length(i1) == 0) {
+    stop("solve: phenotype '", pheno, "' has no samples on one side (degenerate/constant trait) -- cannot compute fold change")
+  }
   m1 <- rowMeans(X[, i1, drop = FALSE], na.rm = TRUE)
   m0 <- rowMeans(X[, i0, drop = FALSE], na.rm = TRUE)
   fc <- m1 - m0
-  rho[is.na(rho)] <- 0
 
-  ## for PHENO nodes 'foldchange' does not make sense
   ii <- grep("PHENO", names(fc))
   if (length(ii)) fc[ii] <- rho[ii]
 
-  ## set node values
+  return(fc)
+
+}
+
+## Attach both node scores to the graph and pick the one that acts as the
+## node value, which is what the edge weighting and the plots read.
+set_node_values <- function(graph, rho, fc, value.type) {
+
   igraph::V(graph)$rho <- rho
   igraph::V(graph)$fc <- fc
-  if (value.type == "rho") {
-    igraph::V(graph)$value <- rho
-  } else {
-    igraph::V(graph)$value <- fc
-  }
+  igraph::V(graph)$value <- if (value.type == "rho") rho else fc
   graph$value.type <- value.type
 
-  ## set edge weights from node values
+  return(graph)
+
+}
+
+## Weight each edge by its correlation, optionally scaled by the geometric
+## mean of the node values at both ends, so that edges between two nodes that
+## respond to the phenotype outweigh equally correlated but flat ones. The
+## SOURCE and SINK edges are bookkeeping and always weigh 1.
+set_edge_weights <- function(graph, fc.weights) {
+
   ww <- 1
   weight.type <- "rho"
   if (fc.weights) {
@@ -91,51 +155,50 @@ solve <- function(obj,
     weight.type <- paste0(weight.type, "*vv")
   }
 
-  ## edge weighting
   ee.rho <- igraph::E(graph)$rho
   ee.rho[is.na(ee.rho)] <- 0.1234
   igraph::E(graph)$weight <- ee.rho * ww
 
-  ## set SINK/SOURCE edges to 1
   if (any(grepl("SINK|SOURCE", igraph::V(graph)$name))) {
     igraph::E(graph)[.to("SINK")]$weight <- 1
     igraph::E(graph)[.from("SOURCE")]$weight <- 1
   }
-
-  if (sp.weight) {
-    sp.wt <- sp_edge_weight(graph, obj$layers)
-    sp.wt <- (sp.wt / max(sp.wt, na.rm = TRUE))^2
-    igraph::E(graph)$weight <- igraph::E(graph)$weight * sp.wt
-    weight.type <- paste0(weight.type, "*sp")
-  }
   graph$weight.type <- weight.type
 
-  ## threshold by minimum weight
-  if (min_rho > 0) {
-    dsel <- which(abs(igraph::E(graph)$weight) < min_rho)
-    igraph::E(graph)$weight[dsel] <- 0
-  }
+  return(graph)
 
-  ## limit edges per connection type
-  if (max_edges > 0) {
-    ewt <- igraph::E(graph)$weight
-    esel <- tapply(
-      seq_along(igraph::E(graph)), igraph::E(graph)$connection_type,
-      function(ii) utils::head(ii[order(-abs(ewt[ii]))], max_edges)
-    )
-    dsel <- setdiff(seq_along(igraph::E(graph)), unlist(esel))
-    igraph::E(graph)$weight[dsel] <- 0
-  }
+}
 
-  ## delete zero edges
-  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
+## Scale the edge weights by how well each edge carries a strong path from
+## SOURCE to SINK, which favours edges of complete cross-layer chains over
+## strong but isolated ones. Needs the helper nodes that add.sink adds.
+add_sp_weight <- function(graph, layers) {
 
-  ## prune vertices
-  if (prune) {
-    ewt <- igraph::E(graph)$weight
-    graph <- igraph::subgraph_from_edges(graph, which(abs(ewt) > 0))
+  if (!all(c("SOURCE", "SINK") %in% igraph::V(graph)$name)) {
+    stop("sp.weight=TRUE requires a model built with add.sink=TRUE")
   }
 
+  sp.wt <- sp_edge_weight(graph, layers)
+  sp.wt <- (sp.wt / max(sp.wt, na.rm = TRUE))^2
+  igraph::E(graph)$weight <- igraph::E(graph)$weight * sp.wt
+  graph$weight.type <- paste0(graph$weight.type, "*sp")
+
+  return(graph)
+
+}
+
+## Zero all but the 'max_edges' strongest edges of each connection type, so
+## that no single type can crowd out the others.
+zero_excess_edges <- function(graph, max_edges) {
+
+  ewt <- igraph::E(graph)$weight
+  esel <- tapply(
+    seq_along(igraph::E(graph)), igraph::E(graph)$connection_type,
+    function(ii) utils::head(ii[order(-abs(ewt[ii]))], max_edges)
+  )
+  dsel <- setdiff(seq_along(igraph::E(graph)), unlist(esel))
+  igraph::E(graph)$weight[dsel] <- 0
+
   return(graph)
 
 }
@@ -205,27 +268,43 @@ multisolve <- function(obj,
                        prune = TRUE,
                        sp.weight = FALSE) {
 
+  if (is.null(obj)) {
+    message("[lasagna::multisolve]: obj is NULL. Please check run of create_model(). Exiting.")
+    return(NULL)
+  }
+
   if (is.null(traits)) traits <- colnames(obj$Y)
   traits <- intersect(traits, colnames(obj$Y))
 
   M <- list()
   V <- list()
   for (ct in traits) {
-    solved <- solve(obj,
-      pheno = ct,
-      min_rho = min_rho,
-      max_edges = max_edges,
-      value.type = value.type,
-      fc.weights = fc.weights,
-      sp.weight = sp.weight,
-      prune = FALSE
+    solved <- tryCatch(
+      solve(obj,
+        pheno = ct,
+        min_rho = min_rho,
+        max_edges = max_edges,
+        value.type = value.type,
+        fc.weights = fc.weights,
+        sp.weight = sp.weight,
+        prune = FALSE
+      ),
+      error = function(e) {
+        warning("multisolve: skipping trait '", ct, "': ", conditionMessage(e), call. = FALSE)
+        NULL
+      }
     )
+    if (is.null(solved)) next
     adj <- igraph::as_adjacency_matrix(solved, attr = "weight")
     M[[ct]] <- adj
     V[[ct]] <- igraph::V(solved)$value
     names(V[[ct]]) <- igraph::V(solved)$name
   }
 
+  if (length(M) == 0) {
+    stop("multisolve: no traits produced a valid solve (all were degenerate/constant)")
+  }
+
   ## RMS adjacency matrix as consensus solution
   M <- lapply(M, function(mat) as.matrix(mat^2))
   avgM <- sqrt(Reduce("+", M) / length(M))
diff --git a/R/utils.R b/R/utils.R
index 70e9b0c..feb1009 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -87,7 +87,7 @@ mofa.merge_data2 <- function(xdata, merge.rows = "prefix", merge.cols = "union")
     A <- xdata[[i]]
     ii <- match(rownames(D), rownames(A))
     jj <- match(colnames(D), colnames(A))
-    A1 <- A[ii, jj]
+    A1 <- A[ii, jj, drop = FALSE]
     nn <- nn + !is.na(A1) * 1
     A1[is.na(A1)] <- 0
     D <- D + A1
diff --git a/man/create_model.Rd b/man/create_model.Rd
index 60437eb..7b5600a 100644
--- a/man/create_model.Rd
+++ b/man/create_model.Rd
@@ -27,6 +27,16 @@ create_model(
 \item{data}{A list with \code{X} (named list of data matrices),
 \code{samples} (data frame), and optionally \code{contrasts}.}
 
+\item{X}{Named list of data matrices, one per layer. Alternative to
+\code{data} when \code{data} is not supplied.}
+
+\item{meta}{Sample phenotype data frame (or contrasts matrix, per
+\code{meta.type}). Alternative to \code{data} when \code{data} is
+not supplied.}
+
+\item{meta.type}{Phenotype type: \code{"pheno"}, \code{"expanded"},
+or \code{"contrasts"}.}
+
 \item{ntop}{Number of top-SD features per layer. Set 0 or NULL
 to keep all.}
 
@@ -47,9 +57,6 @@ to keep all.}
 \item{add.revpheno}{Add reversed phenotype contrasts.}
 
 \item{condition.edges}{Weight edges by phenotype correlation.}
-
-\item{pheno}{Phenotype type: \code{"pheno"}, \code{"expanded"},
-or \code{"contrasts"}.}
 }
 \value{
 A list with components:
@@ -66,3 +73,15 @@ Create LASAGNA multi-layer graph model
 Builds a multi-partite graph from multi-omics data layers.
 Edges weighted by correlation, optionally conditioned on phenotype.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+model$layers
+}
diff --git a/man/lasagna-package.Rd b/man/lasagna-package.Rd
index de9a6c0..8b967f0 100644
--- a/man/lasagna-package.Rd
+++ b/man/lasagna-package.Rd
@@ -6,6 +6,23 @@
 \alias{lasagna-package}
 \title{lasagna: Multi-Layer Graph Analysis for Multi-Omics Data}
 \description{
-LASAGNA (multi-layer graph analysis) builds multi-partite graphs from multi-omics data layers and solves for condition-specific subnetworks. Supports interactive 3D visualization with plotly, visNetwork, and grimon.
+LASAGNA (multi-layer graph analysis) builds multi-partite graphs from multi-omics data layers and solves for condition-specific subnetworks. Edges between layers are weighted by feature-feature correlation and can be conditioned on a phenotype of interest. The package provides interactive 2D and 3D visualisations of the resulting networks using plotly and visNetwork.
+}
+\seealso{
+Useful links:
+\itemize{
+  \item \url{https://github.com/bigomics/lasagna}
+  \item Report bugs at \url{https://github.com/bigomics/lasagna/issues}
+}
+
+}
+\author{
+\strong{Maintainer}: Ivo Kwee \email{ivo.kwee@bigomics.com}
+
+Other contributors:
+\itemize{
+  \item BigOmics Analytics SA [copyright holder, funder]
+}
+
 }
 \keyword{internal}
diff --git a/man/layout_multipartite_3d.Rd b/man/layout_multipartite_3d.Rd
new file mode 100644
index 0000000..fb63dcb
--- /dev/null
+++ b/man/layout_multipartite_3d.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{layout_multipartite_3d}
+\alias{layout_multipartite_3d}
+\title{Compute 3D layout coordinates per layer via dimensionality reduction
+Splits the input data matrix by layer and applies SVD, t-SNE, or UMAP
+to compute 2D (x, y) coordinates within each layer, combined with a
+layer factor as the z-coordinate, for use as a layout by \code{plot_3d}.}
+\usage{
+layout_multipartite_3d(graph, X, clust = c("svd", "tsne", "umap"))
+}
+\arguments{
+\item{graph}{An igraph object (output of \code{solve}).}
+
+\item{X}{Numeric matrix of features (rows, named to match graph vertex
+names) by samples (columns), spanning all layers.}
+
+\item{clust}{Dimensionality reduction method used within each layer:
+\code{"svd"}, \code{"tsne"}, or \code{"umap"}.}
+}
+\value{
+A data frame with columns \code{x}, \code{y} (per-layer 2D
+  layout coordinates) and \code{z} (a factor giving the layer of each
+  vertex), one row per graph vertex, row names matching vertex names.
+}
+\description{
+Compute 3D layout coordinates per layer via dimensionality reduction
+Splits the input data matrix by layer and applies SVD, t-SNE, or UMAP
+to compute 2D (x, y) coordinates within each layer, combined with a
+layer factor as the z-coordinate, for use as a layout by \code{plot_3d}.
+}
diff --git a/man/multisolve.Rd b/man/multisolve.Rd
index aefa818..ec37b86 100644
--- a/man/multisolve.Rd
+++ b/man/multisolve.Rd
@@ -43,3 +43,16 @@ Solve LASAGNA graph for multiple phenotypes
 Solves the graph iteratively for multiple contrasts and returns
 the root-mean-square (RMS) consensus graph.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- multisolve(model, min_rho = 0, prune = FALSE)
+igraph::ecount(g)
+}
diff --git a/man/plot_3d.Rd b/man/plot_3d.Rd
index 72764b6..72f783f 100644
--- a/man/plot_3d.Rd
+++ b/man/plot_3d.Rd
@@ -26,15 +26,39 @@ plot_3d(
 \arguments{
 \item{graph}{An igraph object (output of \code{solve}).}
 
+\item{layout}{Either a matrix/data frame with columns \code{x}, \code{y},
+\code{z} (one row per vertex, row names matching vertex names), or a named
+list of 2-column position matrices per layer. If \code{NULL}, a
+\code{layout} attribute pre-set on \code{graph} (\code{graph$layout <-
+...}) is used before falling back to an automatically computed layout.}
+
+\item{X}{Numeric matrix of features (rows, named to match graph vertex
+names) by samples (columns), used to compute the layout when
+\code{layout} is \code{NULL} or a method name.}
+
 \item{draw_edges}{Logical; draw inter-layer edges.}
 
 \item{num_edges}{Maximum number of edges per layer pair.}
 
 \item{min_rho}{Minimum absolute weight for edges.}
 
+\item{sign_rho}{Which edge signs to draw: \code{"pos"}, \code{"neg"} or
+\code{"both"}.}
+
+\item{cex}{Point size multiplier.}
+
+\item{cex.gamma}{Gamma exponent applied to point sizes.}
+
+\item{edge.cex}{Edge width multiplier.}
+
+\item{color.by}{Vertex colouring: \code{"value"} or \code{"color"}.}
+
+\item{edge_colors}{Two colors for negative and positive edge weights.}
+
 \item{znames}{Named character vector mapping layer codes to display names.}
 
-\item{pos}{Named list of 2-column position matrices per layer.}
+\item{ax}{Integer (mod 4) selecting which plot corner the layer title
+text is placed at.}
 }
 \value{
 A plotly object.
@@ -44,3 +68,22 @@ Plot LASAGNA graph in 3D using plotly
 Wrapper that creates a 3D plotly visualization from a solved
 LASAGNA graph and precomputed 2D positions per layer.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- solve(model, pheno = colnames(model$Y)[1], prune = FALSE)
+## build a simple per-layer 3D layout (random x/y, z = layer)
+vn <- igraph::V(g)$name
+layout <- data.frame(
+  x = runif(length(vn)), y = runif(length(vn)),
+  z = sub(":.*", "", vn),
+  row.names = sub(".*:", "", vn))
+fig <- plot_3d(g, layout = layout, num_edges = 20, min_rho = 0)
+}
diff --git a/man/plot_multipartite.Rd b/man/plot_multipartite.Rd
index c32903b..d071882 100644
--- a/man/plot_multipartite.Rd
+++ b/man/plot_multipartite.Rd
@@ -62,6 +62,10 @@ plot_multipartite(
 
 \item{edge.alpha}{Edge transparency.}
 
+\item{edge.gamma}{Gamma exponent applied to scaled edge widths.}
+
+\item{edge.colors}{Two colors for negative and positive edge weights.}
+
 \item{xdist}{Distance between layers.}
 
 \item{normalize.edges}{Normalize edges per connection type.}
@@ -82,13 +86,32 @@ plot_multipartite(
 
 \item{prune}{Remove disconnected vertices.}
 
+\item{do.plot}{Logical; if \code{FALSE}, compute the layout without drawing.}
+
+\item{color.var}{Vertex colouring: \code{"value"} (by sign) or
+\code{"color"} (use the vertex \code{color} attribute).}
+
 \item{layout}{Layout type: \code{"parallel"} or \code{"hive"}.}
 }
 \value{
-Invisibly returns NULL. Called for side effect (plot).
+Invisibly, a list with the pruned \code{graph} and its
+  \code{layout} coordinate matrix. Called mainly for the plot side effect.
 }
 \description{
 Plot multi-partite graph using base R graphics
 Draws a multi-layer LASAGNA graph as a parallel coordinate-style
 layout using base R \code{igraph::plot}.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- solve(model, pheno = colnames(model$Y)[1], prune = FALSE)
+mp <- plot_multipartite(g, min.rho = 0, ntop = 10)
+}
diff --git a/man/plot_visgraph.Rd b/man/plot_visgraph.Rd
index be624d1..405f707 100644
--- a/man/plot_visgraph.Rd
+++ b/man/plot_visgraph.Rd
@@ -35,6 +35,20 @@ plot_visgraph(
 
 \item{ecex}{Edge width multiplier.}
 
+\item{egamma}{Gamma exponent applied to edge widths.}
+
+\item{color.var}{Vertex colouring: \code{"value"}, \code{"type"}/\code{"layer"},
+or \code{"color"} (use the vertex \code{color} attribute).}
+
+\item{edge.colors}{Two colors for negative and positive edge weights.}
+
+\item{labcex}{Label size multiplier.}
+
+\item{layout}{Optional layout matrix (rows named by vertex). If
+\code{NULL}, a \code{layout} attribute pre-set on \code{graph}
+(\code{graph$layout <- ...}) is used before falling back to the
+default visNetwork layout.}
+
 \item{physics}{Enable physics simulation.}
 }
 \value{
@@ -43,3 +57,16 @@ A visNetwork widget.
 \description{
 Plot LASAGNA graph as interactive network with visNetwork
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- solve(model, pheno = colnames(model$Y)[1], prune = FALSE)
+vis <- plot_visgraph(g, min_rho = 0, ntop = 20)
+}
diff --git a/man/plotlyLasagna.Rd b/man/plotlyLasagna.Rd
index 0bda9b3..1e17cd2 100644
--- a/man/plotlyLasagna.Rd
+++ b/man/plotlyLasagna.Rd
@@ -20,9 +20,16 @@ plotlyLasagna(
 
 \item{znames}{Named character vector for layer display names.}
 
+\item{ax}{Integer (mod 4) selecting which plot corner the layer title
+text is placed at.}
+
 \item{cex}{Point size multiplier.}
 
+\item{edge.cex}{Edge line width multiplier.}
+
 \item{edges}{Optional data frame of edges.}
+
+\item{edge_colors}{Two colors for negative and positive edge weights.}
 }
 \value{
 A plotly object.
@@ -31,3 +38,13 @@ A plotly object.
 Internal plotly builder for LASAGNA 3D plot
 Builds the actual plotly figure from a prepared data frame.
 }
+\examples{
+set.seed(1)
+df <- data.frame(
+  feature = paste0("f", 1:10),
+  x = runif(10), y = runif(10),
+  z = rep(c("gx", "px"), each = 5),
+  color = rnorm(10), size = runif(10),
+  text = paste0("f", 1:10))
+fig <- plotlyLasagna(df)
+}
diff --git a/man/prune_graph.Rd b/man/prune_graph.Rd
index ad36fc3..28e90a3 100644
--- a/man/prune_graph.Rd
+++ b/man/prune_graph.Rd
@@ -52,3 +52,17 @@ Filters vertices and edges by layer, node value, edge weight,
 edge sign, and edge type (inter vs intra). Useful for reducing
 graph complexity before visualization.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- solve(model, pheno = colnames(model$Y)[1], prune = FALSE)
+pruned <- prune_graph(g, ntop = 5, min.rho = 0)
+igraph::vcount(pruned)
+}
diff --git a/man/solve.Rd b/man/solve.Rd
index 6169efa..d821c90 100644
--- a/man/solve.Rd
+++ b/man/solve.Rd
@@ -47,3 +47,16 @@ Given a LASAGNA model and a phenotype (column of Y), computes
 per-node fold change and correlation, then weights edges accordingly.
 The resulting graph is pruned to \code{max_edges} per connection type.
 }
+\examples{
+set.seed(1)
+gx <- matrix(rnorm(20 * 10), 20, 10,
+  dimnames = list(paste0("g", 1:20), paste0("S", 1:10)))
+px <- matrix(rnorm(15 * 10), 15, 10,
+  dimnames = list(paste0("p", 1:15), paste0("S", 1:10)))
+samples <- data.frame(group = rep(c("A", "B"), each = 5),
+  row.names = paste0("S", 1:10))
+data <- list(X = list(gx = gx, px = px), samples = samples)
+model <- create_model(data, ntop = 10, nc = 5)
+g <- solve(model, pheno = colnames(model$Y)[1], max_edges = 50, prune = FALSE)
+igraph::vcount(g)
+}
diff --git a/tests/multipartite-moxbrca.pdf b/tests/multipartite-moxbrca.pdf
deleted file mode 100644
index 1d70947..0000000
Binary files a/tests/multipartite-moxbrca.pdf and /dev/null differ
diff --git a/tests/test-simple.R b/tests/test-simple.R
deleted file mode 100644
index 94f36ef..0000000
--- a/tests/test-simple.R
+++ /dev/null
@@ -1,141 +0,0 @@
-library(data.table)
-library(ggplot2)
-library(igraph)
-library(devtools)
-#load_all("~/Playground/playbase")
-#load_all("~/Projects/WGCNAplus")
-library(WGCNAplus)
-load_all()
-
-## Imports
-gset.rankcor <- playbase::gset.rankcor
-mat2gmt <- playbase::mat2gmt
-ai.ask <- playbase::ai.ask
-ai.create_image_gemini <- playbase::ai.create_image_gemini
-#lasagna.multisolve <- playbase::lasagna.multisolve
-lasagna.multisolve <- multisolve
-
-data <- playbase::mofa.exampledata("brca")
-#data <- playbase::mofa.exampledata("cll")
-
-names(data)
-names(data$X)
-##names(data$X) <- substring(names(data$X),1,2)
-
-ntop=1000;nc=20;add.sink=TRUE;intra=TRUE;use.gmt=FALSE;use.graphite=0;
-fully_connect=FALSE;add.revpheno=TRUE;condition.edges=1
-
-obj <- lasagna::create_model(data, meta.type="pheno", ntop=1000, nc=10,
-  add.sink=TRUE, intra=FALSE, fully_connect=FALSE, add.revpheno=TRUE,
-  condition.edges=1)
-names(obj)
-
-## color by WGCNA clustering
-if(0) {
-
-  names(data$X)
-  wgcna <- computeWGCNA_multiomics(data$X, data$samples,
-    power=12, minmodsize = 3, minKME=0.1, mergeCutHeight = 0.5)
-  names(wgcna)
-  names(wgcna$layers)
-   
-  par(mfrow=c(1,1),cex=1)
-  WGCNAplus::plotMultiDendroAndColors(
-    wgcna, marAll=c(2,7,3,1), #multi=TRUE, 
-    show.traits=1, show.contrasts=1,
-    show.kme=0, use.tree=0,
-    colorHeight = 0.5, main=""
-  )
-  
-  wgcna$me.genes
-  table(wgcna$me.colors)
-  head(wgcna$me.colors)
-  
-  V(obj$graph)$color <- wgcna$me.colors[V(obj$graph)$name]
-  ##V(obj$graph)$color <- substring(wgcna$me.colors[V(obj$graph)$name],3,99)
-  V(obj$graph)$color[is.na(V(obj$graph)$color)] <- "red"
-  table(V(obj$graph)$color)
-}
-
-wgcna$me.genes
-
-## solve the graph for a certain phenotype
-colnames(obj$Y)
-pheno = "activated=act"
-pheno = "condition=Her2"
-pheno = "condition=LumA"
-graph <- lasagna::solve(obj, pheno, min_rho=0.01, max_edges=1000,
-  value="rho", sp.weight=1, prune=FALSE) 
-graph
-
-## graph <- lasagna::multisolve(obj, pheno, min_rho=0.01, max_edges=1000,
-##   value="rho", sp.weight=1, prune=FALSE) 
-## graph
-
-## prune graph for cleaner plotting
-pdf("multipartite-moxbrca.pdf", w=14, h=8)
-par(mfrow=c(1,1), mar=c(1,1,1,1)*0)
-mp <- lasagna::plot_multipartite(
-  graph,
-  min.rho = 0.3,
-  ntop = 50,
-  xdist = 1,
-  color.var = "color",
-  labpos = c(2,2,4,4),
-  cex.label = 0.8,
-  vx.cex = 1.1,
-  edge.cex = 1.5,
-  edge.alpha = 0.4,
-  edge.sign = "both",
-  edge.type = "inter",
-  edge.gamma = 2,
-  yheight = 0.99,
-  normalize.edges = 1,
-  strip.prefix = TRUE,
-  prune = 0
-) 
-dev.off()
-
-names(mp)
-mp$graph
-table(V(mp$graph)$color)
-
-## interactive multipartite
-M <- mp$layout
-vis <- plot_visgraph(mp$graph, layers=NULL, ntop=-1,
-  min_rho=0.2, ecex=3, vcex=3, labcex=1, egamma=2,
-  color.var="color", mst=0, layout=M, physics=0) 
-vis
-
-## MST - FR layout (NEED RETHINK)
-require(visNetwork)
-mst <- mp$graph
-##ew <- 1 / pmax(E(mp$graph)$weight,0)
-ew <- 1 / (1e-8 + abs(E(mp$graph)$weight)**2)
-mst <- igraph::mst(mp$graph, weights=ew)
-vis <- plot_visgraph(mst, layers=NULL, ntop=-1,
-  min_rho=0.1, ecex=3, egamma=2, vcex=2, labcex=1, 
-  mst=1, color.var="color", layout=NULL, physics=TRUE) 
-vis
-
-## hierarchical layout
-vis %>% visHierarchicalLayout()
-
-##-------------------------------------------------
-## 3D lasagna
-##-------------------------------------------------
-##source("~/Playground/playbase/dev/include.R", chdir=TRUE)
-load_all("..")
-
-#xpos <- layout_multipartite_3d(graph, obj$X, clust='svd')
-xpos <- layout_multipartite_3d(graph, obj$X, clust='tsne')
-#xpos <- layout_multipartite_3d(graph, obj$X, clust='umap')
-head(xpos)
-
-plot_3d(graph, layout=xpos, draw_edges=TRUE,
-  color.by="color", min_rho=0.0, sign_rho="both",
-  cex=0.9, cex.gamma=0.5, num_edges=200, znames=NULL) 
-  
-##-------------------------------------------------
-##-------------------------------------------------
-##-------------------------------------------------
diff --git a/tests/testthat/helper-toy.R b/tests/testthat/helper-toy.R
new file mode 100644
index 0000000..5cfddc1
--- /dev/null
+++ b/tests/testthat/helper-toy.R
@@ -0,0 +1,34 @@
+## Shared fixtures for the plotting/graph/solve tests: small toy models
+## built once at suite load and reused across test_that() blocks so the
+## suite stays fast.
+
+make_toy_data <- function(n = 10, p1 = 20, p2 = 15, seed = 42) {
+  set.seed(seed)
+  gx <- matrix(rnorm(p1 * n), nrow = p1, ncol = n)
+  rownames(gx) <- paste0("gene", 1:p1)
+  colnames(gx) <- paste0("S", 1:n)
+
+  px <- matrix(rnorm(p2 * n), nrow = p2, ncol = n)
+  rownames(px) <- paste0("prot", 1:p2)
+  colnames(px) <- paste0("S", 1:n)
+
+  samples <- data.frame(
+    group = factor(rep(c("A", "B"), each = n / 2)),
+    row.names = paste0("S", 1:n)
+  )
+
+  list(X = list(gx = gx, px = px), samples = samples)
+}
+
+toy_model <- suppressMessages(
+  create_model(make_toy_data(), meta.type = "pheno", ntop = 10, nc = 5)
+)
+toy_pheno <- colnames(toy_model$Y)[1]
+toy_graph <- solve(toy_model, pheno = toy_pheno, max_edges = 50, prune = FALSE)
+
+## a layer reduced to a single feature, the shape that used to crash
+## layout_multipartite_3d's per-layer svd
+single_feature_model <- suppressMessages(
+  create_model(make_toy_data(n = 8, p1 = 20, p2 = 1, seed = 7),
+    meta.type = "pheno", ntop = 10, nc = 5)
+)
diff --git a/tests/testthat/test-graph.R b/tests/testthat/test-graph.R
new file mode 100644
index 0000000..5b5fb6f
--- /dev/null
+++ b/tests/testthat/test-graph.R
@@ -0,0 +1,30 @@
+test_that("prune_graph filters by layers and ntop", {
+  skip_on_cran()
+  pruned <- prune_graph(toy_graph, layers = "gx", ntop = 5, min.rho = 0)
+  expect_true(all(igraph::V(pruned)$layer == "gx"))
+  expect_lte(igraph::vcount(pruned), 5)
+})
+
+test_that("prune_graph filters edges by sign", {
+  skip_on_cran()
+  pruned_pos <- prune_graph(toy_graph, ntop = 10, min.rho = 0, edge.sign = "pos", prune = FALSE)
+  expect_true(all(igraph::E(pruned_pos)$weight > 0))
+
+  pruned_neg <- prune_graph(toy_graph, ntop = 10, min.rho = 0, edge.sign = "neg", prune = FALSE)
+  expect_true(all(igraph::E(pruned_neg)$weight < 0))
+})
+
+test_that("prune_graph filters edges by type", {
+  skip_on_cran()
+  pruned_inter <- prune_graph(toy_graph, ntop = 10, min.rho = 0, edge.type = "inter", prune = FALSE)
+  expect_true(all(grepl("->", igraph::E(pruned_inter)$connection_type)))
+
+  pruned_intra <- prune_graph(toy_graph, ntop = 10, min.rho = 0, edge.type = "intra", prune = FALSE)
+  expect_true(all(!grepl("->", igraph::E(pruned_intra)$connection_type)))
+})
+
+test_that("prune_graph applies a minimum rho threshold", {
+  skip_on_cran()
+  pruned <- prune_graph(toy_graph, ntop = 10, min.rho = 0.3, prune = FALSE)
+  expect_true(all(abs(igraph::E(pruned)$weight) >= 0.3))
+})
diff --git a/tests/testthat/test-model.R b/tests/testthat/test-model.R
index d791830..51e0bce 100644
--- a/tests/testthat/test-model.R
+++ b/tests/testthat/test-model.R
@@ -24,7 +24,7 @@ test_that("create_model works with minimal data", {
     samples = samples
   )
 
-  model <- create_model(data, pheno = "pheno", ntop = 10, nc = 5)
+  model <- create_model(data, meta.type = "pheno", ntop = 10, nc = 5)
 
   expect_true(is.list(model))
   expect_true("graph" %in% names(model))
@@ -59,7 +59,7 @@ test_that("solve works on created model", {
     samples = samples
   )
 
-  model <- create_model(data, pheno = "pheno", ntop = 10, nc = 5)
+  model <- create_model(data, meta.type = "pheno", ntop = 10, nc = 5)
   pheno <- colnames(model$Y)[1]
 
   solved <- solve(model, pheno = pheno, max_edges = 50)
diff --git a/tests/testthat/test-plotting.R b/tests/testthat/test-plotting.R
new file mode 100644
index 0000000..67cf701
--- /dev/null
+++ b/tests/testthat/test-plotting.R
@@ -0,0 +1,123 @@
+## layout_multipartite_3d ---------------------------------------------------
+
+test_that("layout_multipartite_3d default clust argument does not crash", {
+  skip_on_cran()
+  out <- layout_multipartite_3d(toy_model$graph, toy_model$X)
+  expect_s3_class(out, "data.frame")
+  expect_true(all(c("x", "y", "z") %in% colnames(out)))
+  expect_equal(nrow(out), igraph::vcount(toy_model$graph))
+})
+
+test_that("layout_multipartite_3d runs with an explicit clust value", {
+  skip_on_cran()
+  out <- layout_multipartite_3d(toy_model$graph, toy_model$X, clust = "svd")
+  expect_s3_class(out, "data.frame")
+  expect_equal(rownames(out), igraph::V(toy_model$graph)$name)
+})
+
+test_that("layout_multipartite_3d handles a single-feature layer", {
+  skip_on_cran()
+  out <- layout_multipartite_3d(
+    single_feature_model$graph, single_feature_model$X, clust = "svd"
+  )
+  expect_s3_class(out, "data.frame")
+  px_row <- grep("^px:", rownames(out))
+  expect_length(px_row, 1)
+  expect_equal(unname(unlist(out[px_row, c("x", "y")])), c(0, 0))
+})
+
+## plot_multipartite ---------------------------------------------------------
+
+test_that("plot_multipartite runs and returns the pruned graph and layout", {
+  skip_on_cran()
+  out <- plot_multipartite(toy_graph, min.rho = 0, ntop = 10, do.plot = FALSE)
+  expect_true(is.list(out))
+  expect_true(igraph::is_igraph(out$graph))
+  expect_true(is.matrix(out$layout))
+})
+
+test_that("plot_multipartite errors clearly when the layers filter matches nothing", {
+  skip_on_cran()
+  expect_error(
+    plot_multipartite(toy_graph, layers = "nonexistent", do.plot = TRUE),
+    "no nodes/edges remain"
+  )
+  expect_error(
+    plot_multipartite(toy_graph, layers = "nonexistent", do.plot = FALSE),
+    "no nodes/edges remain"
+  )
+})
+
+## plot_3d --------------------------------------------------------------------
+
+test_that("plot_3d runs with clust = 'svd'", {
+  skip_on_cran()
+  fig <- plot_3d(toy_graph, layout = "svd", X = toy_model$X, num_edges = 5, min_rho = 0)
+  expect_s3_class(fig, "plotly")
+})
+
+test_that("plot_3d uses a layout pre-attached to the graph", {
+  skip_on_cran()
+  g <- toy_graph
+  g$layout <- layout_multipartite_3d(g, toy_model$X, clust = "svd")
+  fig <- plot_3d(g, num_edges = 5, min_rho = 0)
+  expect_s3_class(fig, "plotly")
+})
+
+test_that("plot_3d errors clearly when neither X nor layout is given", {
+  skip_on_cran()
+  expect_error(plot_3d(toy_graph), "must provide X or layout")
+})
+
+test_that("plot_3d errors clearly on an incomplete layout", {
+  skip_on_cran()
+  incomplete_X <- toy_model$X[-1, , drop = FALSE]
+  expect_error(
+    plot_3d(toy_graph, layout = "svd", X = incomplete_X),
+    "incomplete layout"
+  )
+})
+
+test_that("plot_3d errors clearly on an invalid layout string", {
+  skip_on_cran()
+  expect_error(
+    plot_3d(toy_graph, layout = "badmethod", X = toy_model$X),
+    "invalid layout"
+  )
+})
+
+## plotlyLasagna ---------------------------------------------------------------
+
+test_that("plotlyLasagna renders a manually constructed data frame", {
+  skip_on_cran()
+  set.seed(1)
+  df <- data.frame(
+    feature = paste0("f", 1:10),
+    x = runif(10), y = runif(10),
+    z = rep(c("gx", "px"), each = 5),
+    color = rnorm(10), size = runif(10),
+    text = paste0("f", 1:10)
+  )
+  fig <- plotlyLasagna(df)
+  expect_s3_class(fig, "plotly")
+})
+
+## plot_visgraph -----------------------------------------------------------
+
+test_that("plot_visgraph returns a visNetwork widget", {
+  skip_on_cran()
+  vis <- plot_visgraph(toy_graph, min_rho = 0, ntop = 20)
+  expect_s3_class(vis, "visNetwork")
+  expect_s3_class(vis, "htmlwidget")
+})
+
+test_that("plot_visgraph supports a non-default color.var and an explicit layout", {
+  skip_on_cran()
+  pos <- cbind(
+    x = stats::rnorm(igraph::vcount(toy_graph)),
+    y = stats::rnorm(igraph::vcount(toy_graph))
+  )
+  rownames(pos) <- igraph::V(toy_graph)$name
+  vis <- plot_visgraph(toy_graph, min_rho = 0, ntop = 20, color.var = "layer", layout = pos)
+  expect_s3_class(vis, "visNetwork")
+})
diff --git a/tests/testthat/test-solve2.R b/tests/testthat/test-solve2.R
new file mode 100644
index 0000000..3a837fc
--- /dev/null
+++ b/tests/testthat/test-solve2.R
@@ -0,0 +1,7 @@
+test_that("multisolve runs across all traits and returns a valid igraph", {
+  skip_on_cran()
+  g <- multisolve(toy_model, min_rho = 0, prune = FALSE)
+  expect_true(igraph::is_igraph(g))
+  expect_equal(igraph::vcount(g), igraph::vcount(toy_model$graph))
+  expect_true("value" %in% igraph::vertex_attr_names(g))
+})