Optimize

.Rbuildignore

@@ -33,3 +33,5 @@ vignettes/.RData
 ^CODE_OF_CONDUCT\.md$
 ^_pkgdown\.yml$
 ^codemeta\.json$
+^\.positai$
+^\.claude$

.github/workflows/memcheck.yml

@@ -1,5 +1,3 @@
-# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
-# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
 on:
   push:
     branches:
@@ -24,7 +22,7 @@ name: mem-check
 
 jobs:
   mem-check:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-24.04
 
     name: valgrind ${{ matrix.config.test }}
 
@@ -36,67 +34,11 @@ jobs:
           - {test: 'examples'}
 
     env:
-      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
       _R_CHECK_FORCE_SUGGESTS_: false
-      RSPM: https://packagemanager.rstudio.com/cran/__linux__/noble/latest
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
 
     steps:
-      - uses: actions/checkout@v5
+      - uses: ms609/actions/memcheck@main
         with:
+          test: ${{ matrix.config.test }}
           submodules: true
-
-      - uses: r-lib/actions/setup-r@v2
-        with:
-          r-version: release # CRAN uses devel, but takes ages to load deps.
-
-      - name: Install apt packages
-        run: |
-          sudo apt-get install valgrind texlive-latex-base
-
-      - uses: r-lib/actions/setup-pandoc@v2
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = c("soft",
-              "Config/Needs/github-actions", "Config/Needs/memcheck")),
-            ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Cache R packages
-        uses: actions/cache@v4
-        with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
-      - name: Install system dependencies
-        run: |
-          while read -r cmd
-          do
-            eval sudo $cmd
-          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
-
-      - name: Install dependencies
-        run: |
-          remotes::install_github('emmanuelparadis/ape')
-          if ("${{ matrix.config.test }}" == 'tests') {
-            remotes::install_github('r-lib/brio#20') # brio mem leak
-          }
-          remotes::install_deps(dependencies = c("soft", "Config/Needs/github-actions", "Config/Needs/memcheck"))
-        shell: Rscript {0}
-
-      - name: Install this package
-        run: |
-          cd ..
-          R CMD build --no-manual --no-resave-data Rogue
-          R CMD INSTALL Rogue*.tar.gz
-          cd Rogue
-
-      - name: valgrind - memcheck ${{ matrix.config.test }}
-        run: |
-          R -d "valgrind --tool=memcheck --leak-check=full \
-          --errors-for-leak-kinds=definite --error-exitcode=1" \
-          --vanilla < memcheck/${{ matrix.config.test }}.R

R/RogueTaxa.R

@@ -159,13 +159,17 @@ RogueTaxa <- function(trees,
                                mreOptimization = mreOptimization,
                                threshold = threshold, verbose = verbose),
                    Roguehalla(trees, dropsetSize = dropsetSize, info = "phylo",
-                              p = p, neverDrop = neverDrop),
+                              p = p, neverDrop = neverDrop,
+                              .prepared = TRUE),
                    Roguehalla(trees, dropsetSize = dropsetSize, info = "clust",
-                              p = p, neverDrop = neverDrop),
+                              p = p, neverDrop = neverDrop,
+                              .prepared = TRUE),
                    QuickRogue(trees, info = "phylo", p = p,
-                              neverDrop = neverDrop),
+                              neverDrop = neverDrop,
+                              .prepared = TRUE),
                    QuickRogue(trees, info = "clust", p = p,
-                              neverDrop = neverDrop)
+                              neverDrop = neverDrop,
+                              .prepared = TRUE)
   )
 
   # Format return value

R/SPIC.R

@@ -12,6 +12,8 @@
 #' @examples
 #'
 #' QuickRogue(trees, fullSeq = TRUE)
+#' @param .prepared Logical; if `TRUE`, skip internal tree preparation.
+#' For internal use.
 #' @importFrom cli cli_progress_bar cli_progress_update cli_progress_done
 #' @importFrom TreeDist ConsensusInfo SplitwiseInfo ClusteringInfo
 #' @importFrom fastmatch %fin%
@@ -22,7 +24,8 @@ QuickRogue <- function(trees,
                        p = 0.5,
                        log = TRUE, average = "median", deviation = "mad",
                        neverDrop, fullSeq = FALSE,
-                       parallel = FALSE) {
+                       parallel = FALSE,
+                       .prepared = FALSE) {
   if (!is.na(pmatch(tolower(info), "spic"))) {
     info <- "phylogenetic"
   } else if (!is.na(pmatch(tolower(info), "scic"))) {
@@ -41,7 +44,7 @@ QuickRogue <- function(trees,
       stop("`trees` must be a list of `phylo` objects")
     }
   }
-  trees <- .PrepareTrees(trees)
+  if (!.prepared) trees <- .PrepareTrees(trees)
   nTip <- NTip(trees[[1]])
 
   tr <- trees
@@ -54,10 +57,17 @@ QuickRogue <- function(trees,
   TotalInfo <- switch(pmatch(info, c("phylogenetic", "clustering")),
                       SplitwiseInfo, ClusteringInfo)
 
+  # Precompute upper bounds: max info for a fully-resolved tree with n leaves
+  bestPossible <- vapply(
+    nTip - 1L - seq_len(nDrops),
+    function(n) TotalInfo(PectinateTree(n)),
+    double(1)
+  )
+
   cli_progress_bar("Drop leaf", total = nDrops * (nDrops + 1L) / 2,
                    .auto_close = FALSE)
   for (i in 1L + seq_len(nDrops)) {
-    bestPossibleNext <- TotalInfo(PectinateTree(nTip - i))
+    bestPossibleNext <- bestPossible[i - 1L]
     bestYet <- max(score, na.rm = TRUE)
     if (bestPossibleNext < bestYet) {
       # message("Broken out: can't attain ", signif(score[i]), " bits with ",
@@ -143,7 +153,7 @@ QuickRogue <- function(trees,
 #' @importFrom TreeTools DropTipPhylo SplitFrequency Preorder RenumberTips
 #' @importFrom utils combn
 Roguehalla <- function(trees, dropsetSize = 1, info = "phylogenetic",
-                        p = 0.5, neverDrop) {
+                        p = 0.5, neverDrop, .prepared = FALSE) {
   if (!inherits(trees, "multiPhylo")) {
     if (inherits(trees, "phylo")) {
       return(data.frame(num = 0,
@@ -158,8 +168,10 @@ Roguehalla <- function(trees, dropsetSize = 1, info = "phylogenetic",
     }
     trees <- structure(trees, class = "multiPhylo")
   }
-  trees <- lapply(trees, RenumberTips, trees[[1]])
-  trees <- lapply(trees, Preorder)
+  if (!.prepared) {
+    trees <- lapply(trees, RenumberTips, trees[[1]])
+    trees <- lapply(trees, Preorder)
+  }
   startTrees <- trees
   labels <- startTrees[[1]]$tip.label
 

R/stability.R

@@ -133,31 +133,64 @@ TipInstability <- function(trees, log = TRUE, average = "mean",
     nTip <- NTip(trees[[1]])
   }
 
-  dists <- vapply(trees, GraphGeodesic, double(nTip * nTip),
-                  nTip = nTip, log = log, asMatrix = FALSE)
+  lt_idx <- which(lower.tri(matrix(TRUE, nTip, nTip)))
+
+  nEdge <- nrow(trees[[1]]$edge)
+  nNode <- trees[[1]]$Nnode
+  # Batch C path requires uniform tree dimensions and log = TRUE
+  useBatch <- isTRUE(log) &&
+    all(vapply(trees, function(tr) nrow(tr$edge) == nEdge, logical(1)))
+
+  if (useBatch) {
+    # Batch C call: returns lower-triangle entries directly (nPairs × nTree)
+    # Ensure preorder (the per-tree GraphGeodesic path does this internally)
+    trees <- lapply(trees, Preorder)
+    parent_all <- integer(nEdge * length(trees))
+    child_all <- integer(nEdge * length(trees))
+    for (k in seq_along(trees)) {
+      rng <- (k - 1L) * nEdge + seq_len(nEdge)
+      parent_all[rng] <- trees[[k]]$edge[, 1] - 1L
+      child_all[rng] <- trees[[k]]$edge[, 2] - 1L
+    }
+    dists_lt <- matrix(
+      .Call(`LOG_GRAPH_GEODESIC_MULTI`,
+            n_tip = as.integer(nTip),
+            n_node = as.integer(nNode),
+            parent = as.integer(parent_all),
+            child = as.integer(child_all),
+            n_edge = as.integer(nEdge),
+            n_tree = as.integer(length(trees))),
+      ncol = length(trees)
+    )
+  } else {
+    dists <- vapply(trees, GraphGeodesic, double(nTip * nTip),
+                    nTip = nTip, log = log, asMatrix = FALSE)
+    dists_lt <- dists[lt_idx, , drop = FALSE]
+  }
 
   whichDev <- pmatch(tolower(deviation), c("sd", "mad"))
   if (is.na(whichDev)) {
     stop("`deviation` must be 'sd' or 'mad'")
   }
-  devs <- matrix(switch(whichDev,
-                        rowVars(dists, std = TRUE, parallel = parallel),
-                        rowMads(dists, parallel = parallel)),
-                 nTip, nTip)
-  devs[is.nan(devs)] <- 0 # rowVars returns NaN instead of 0
-  #diag(devs) <- 0 # Faster than setting to NA, then using rowMeans(rm.na = TRUE)
-
+  devs_lt <- switch(whichDev,
+                    rowVars(dists_lt, std = TRUE, parallel = parallel),
+                    rowMads(dists_lt, parallel = parallel))
+  devs_lt[is.nan(devs_lt)] <- 0
 
   whichAve <- pmatch(tolower(average), c("mean", "median"))
   if (is.na(whichAve)) {
     stop("`average` must be 'mean' or 'median'")
   }
-  aves <- matrix(switch(whichAve, rowmeans, rowMedians)(dists), nTip, nTip)
+  aves_lt <- switch(whichAve, rowmeans, rowMedians)(dists_lt)
+  meanAve <- mean(aves_lt)
 
-  relDevs <- devs / mean(aves[lower.tri(aves)])
+  # Reconstruct symmetric deviation matrix from lower triangle
+  devs <- matrix(0, nTip, nTip)
+  devs[lt_idx] <- devs_lt
+  devs <- devs + t(devs)
 
   setNames(
-    Rfast::rowmeans(relDevs), # Faster than Rfast::colmeans
+    Rfast::rowmeans(devs) / meanAve, # faster than Rfast::colmeans
     TipLabels(trees[[1]])
   )
 }

src/Rogue_init.c

@@ -9,10 +9,12 @@ extern SEXP RogueNaRok(SEXP, SEXP, SEXP, SEXP, SEXP,
                        SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP GRAPH_GEODESIC(SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP LOG_GRAPH_GEODESIC(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP LOG_GRAPH_GEODESIC_MULTI(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 
 static const R_CallMethodDef callMethods[] = {
   {"GRAPH_GEODESIC", (DL_FUNC) &GRAPH_GEODESIC, 5},
   {"LOG_GRAPH_GEODESIC", (DL_FUNC) &LOG_GRAPH_GEODESIC, 5},
+  {"LOG_GRAPH_GEODESIC_MULTI", (DL_FUNC) &LOG_GRAPH_GEODESIC_MULTI, 6},
   {"RogueNaRok", (DL_FUNC) &RogueNaRok, 10},
   {NULL, NULL, 0}
 };

src/graph_geodesic.c

@@ -111,3 +111,49 @@ SEXP LOG_GRAPH_GEODESIC(SEXP n_tip, SEXP n_node, SEXP parent, SEXP child,
   UNPROTECT(2);
   return(RESULT);
 }
+
+// Batch version: compute log-geodesics for multiple trees in one call.
+// Returns lower-triangle entries only (nPairs x nTree matrix).
+// Reuses a single interim buffer across all trees.
+//
+// parent_all, child_all: concatenated edge arrays for all trees
+//   (each tree contributes n_edge consecutive entries)
+// n_tree: number of trees
+SEXP LOG_GRAPH_GEODESIC_MULTI(SEXP n_tip, SEXP n_node, SEXP parent_all,
+                              SEXP child_all, SEXP n_edge, SEXP n_tree) {
+  const int
+    n_tips = INTEGER(n_tip)[0],
+    n_nodes = INTEGER(n_node)[0],
+    all_nodes = n_tips + n_nodes,
+    n_edges = INTEGER(n_edge)[0],
+    n_trees = INTEGER(n_tree)[0],
+    n_pairs = n_tips * (n_tips - 1) / 2
+  ;
+
+  SEXP RESULT = PROTECT(allocVector(REALSXP, (R_xlen_t)n_pairs * n_trees));
+  SEXP INTERIM = PROTECT(allocVector(INTSXP, all_nodes * all_nodes));
+  double *result = REAL(RESULT);
+  int *interim = INTEGER(INTERIM);
+  const int *par_all = INTEGER(parent_all);
+  const int *ch_all = INTEGER(child_all);
+
+  for (int t = 0; t < n_trees; ++t) {
+    const int *par = par_all + (R_xlen_t)t * n_edges;
+    const int *ch = ch_all + (R_xlen_t)t * n_edges;
+
+    graph_geodesic_phylo(&n_tips, &n_nodes, par, ch, &n_edges,
+                         &all_nodes, interim);
+
+    // Extract lower triangle (row > col in column-major order)
+    double *res_col = result + (R_xlen_t)t * n_pairs;
+    int pair_idx = 0;
+    for (int j = 0; j < n_tips - 1; ++j) {
+      for (int i = j + 1; i < n_tips; ++i) {
+        res_col[pair_idx++] = lg[interim[j + all_nodes * i]];
+      }
+    }
+  }
+
+  UNPROTECT(2);
+  return(RESULT);
+}

tests/testthat/test-stability.R

@@ -14,6 +14,17 @@ test_that("TipInstability() null output", {
                setNames(rep(0, 4), paste0("t", 1:4)))
 })
 
+test_that("TipInstability() batch path handles double edge matrices", {
+  trees <- as.phylo(0:3, 6)
+  # Coerce edge matrices to double (as phylo objects sometimes store them)
+  trees <- lapply(trees, function(tr) {
+    storage.mode(tr$edge) <- "double"
+    tr
+  })
+  class(trees) <- "multiPhylo"
+  expect_no_error(TipInstability(trees, log = TRUE))
+})
+
 test_that("GraphGeodesic() works", {
   Test <- function(tr) {
     tr <- Preorder(tr)