Fix LAPJV alignment regression + coverage for mutual_clustering_impl.h

lap_impl.h: add __attribute__((optimize("align-functions=64",
"align-loops=16"))) on lap() (GCC only) to stabilise instruction
alignment across TU layout changes.

cost_matrix.h: add setWithTranspose() and markTransposed() methods.

lap.cpp: use combined fill+transpose in lapjv() to eliminate the extra
makeTranspose() O(n²) pass introduced by the header refactor.

tree_distance_functions.cpp: add cpp_mci_impl_score() wrapper that calls
TreeDist::mutual_clustering_score() from the installable header, covering
find_exact_matches_raw() and the MCI score computation.

test-mci_impl.R: exercises exact-match early exit and partial-match + LAP
paths through the header implementation.

Author: ms609

R/RcppExports.R

@@ -267,6 +267,10 @@ cpp_transfer_dist_cross_pairs <- function(splits_a, splits_b, n_tip, scale, n_th
     .Call(`_TreeDist_cpp_transfer_dist_cross_pairs`, splits_a, splits_b, n_tip, scale, n_threads)
 }
 
+cpp_mci_impl_score <- function(x, y, n_tips) {
+    .Call(`_TreeDist_cpp_mci_impl_score`, x, y, n_tips)
+}
+
 cpp_robinson_foulds_distance <- function(x, y, nTip) {
     .Call(`_TreeDist_cpp_robinson_foulds_distance`, x, y, nTip)
 }

inst/include/TreeDist/cost_matrix.h

@@ -90,6 +90,16 @@ class CostMatrix {
     return &t_data_[static_cast<std::size_t>(i) * dim8_];
   }
 
+  // Write a value to both the data and transposed-data arrays.
+  // After filling the entire matrix this way, call markTransposed()
+  // to skip the lazy transpose in makeTranspose().
+  void setWithTranspose(lap_row i, lap_col j, cost value) {
+    data_[static_cast<std::size_t>(i) * dim8_ + j] = value;
+    t_data_[static_cast<std::size_t>(j) * dim8_ + i] = value;
+  }
+
+  void markTransposed() noexcept { transposed_ = true; }
+
   // ---- Transpose ----
 
   void makeTranspose() noexcept {

inst/include/TreeDist/lap_impl.h

@@ -32,6 +32,13 @@ namespace detail {
   }
 } // namespace detail
 
+// The LAP hot loops are sensitive to instruction alignment — even unrelated
+// changes to other code in the same translation unit can shift layout enough
+// to cause measurable regressions.  Force 64-byte function alignment and
+// 16-byte loop alignment to stabilise performance across builds.
+#if defined(__GNUC__) && !defined(__clang__)
+__attribute__((optimize("align-functions=64", "align-loops=16")))
+#endif
 cost lap(const lap_row dim,
          CostMatrix& input_cost,
          std::vector<lap_col>& rowsol,

src/RcppExports.cpp

@@ -611,6 +611,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cpp_mci_impl_score
+double cpp_mci_impl_score(const Rcpp::RawMatrix& x, const Rcpp::RawMatrix& y, int n_tips);
+RcppExport SEXP _TreeDist_cpp_mci_impl_score(SEXP xSEXP, SEXP ySEXP, SEXP n_tipsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const Rcpp::RawMatrix& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::RawMatrix& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< int >::type n_tips(n_tipsSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_mci_impl_score(x, y, n_tips));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cpp_robinson_foulds_distance
 List cpp_robinson_foulds_distance(const RawMatrix& x, const RawMatrix& y, const IntegerVector& nTip);
 RcppExport SEXP _TreeDist_cpp_robinson_foulds_distance(SEXP xSEXP, SEXP ySEXP, SEXP nTipSEXP) {
@@ -753,6 +766,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_TreeDist_cpp_transfer_dist_scored", (DL_FUNC) &_TreeDist_cpp_transfer_dist_scored, 4},
     {"_TreeDist_cpp_transfer_dist_all_pairs", (DL_FUNC) &_TreeDist_cpp_transfer_dist_all_pairs, 4},
     {"_TreeDist_cpp_transfer_dist_cross_pairs", (DL_FUNC) &_TreeDist_cpp_transfer_dist_cross_pairs, 5},
+    {"_TreeDist_cpp_mci_impl_score", (DL_FUNC) &_TreeDist_cpp_mci_impl_score, 3},
     {"_TreeDist_cpp_robinson_foulds_distance", (DL_FUNC) &_TreeDist_cpp_robinson_foulds_distance, 3},
     {"_TreeDist_cpp_robinson_foulds_info", (DL_FUNC) &_TreeDist_cpp_robinson_foulds_info, 3},
     {"_TreeDist_cpp_matching_split_distance", (DL_FUNC) &_TreeDist_cpp_matching_split_distance, 3},

src/lap.cpp

@@ -50,16 +50,25 @@ Rcpp::List lapjv(Rcpp::NumericMatrix &x, Rcpp::NumericVector &maxX) {
   std::vector<lap_row> colsol(max_dim);
 
   // Build cost matrix from R's column-major NumericMatrix.
+  // Fill data and transpose simultaneously to avoid a separate
+  // makeTranspose() pass in lap().
   cost_matrix input(max_dim);
   const double* __restrict__ src_data = REAL(x);
   for (lap_row r = 0; r < n_row; ++r) {
     for (lap_col c = 0; c < n_col; ++c) {
-      input(r, c) = static_cast<cost>(
-        src_data[static_cast<std::size_t>(c) * n_row + r] * scale_factor);
+      input.setWithTranspose(r, c, static_cast<cost>(
+        src_data[static_cast<std::size_t>(c) * n_row + r] * scale_factor));
+    }
+    for (lap_col c = n_col; c < max_dim; ++c) {
+      input.setWithTranspose(r, c, max_score);
+    }
+  }
+  for (lap_row r = n_row; r < max_dim; ++r) {
+    for (lap_col c = 0; c < max_dim; ++c) {
+      input.setWithTranspose(r, c, max_score);
     }
-    input.padRowAfterCol(r, n_col, max_score);
   }
-  input.padAfterRow(n_row, max_score);
+  input.markTransposed();
 
   cost score = lap(max_dim, input, rowsol, colsol);
 

src/tree_distance_functions.cpp

@@ -1,3 +1,4 @@
+#include <TreeTools/SplitList.h>
 #include <Rcpp/Lightest>
 
 // Provide the MCI table definitions and implementation in this TU.
@@ -11,3 +12,36 @@ __attribute__((constructor))
   void initialize_ldf() {
     TreeDist::init_lg2_tables(SL_MAX_TIPS);
   }
+
+// Thin wrapper that exercises the installable-header version of
+// mutual_clustering_score(), for test coverage and downstream validation.
+// [[Rcpp::export]]
+double cpp_mci_impl_score(const Rcpp::RawMatrix& x,
+                          const Rcpp::RawMatrix& y,
+                          int n_tips) {
+  using TreeTools::SplitList;
+
+  const SplitList a(x);
+  const SplitList b(y);
+  TreeDist::LapScratch scratch;
+
+  // Build arrays matching the header's raw-pointer API types.
+  std::vector<const splitbit*> a_ptrs(a.n_splits);
+  std::vector<const splitbit*> b_ptrs(b.n_splits);
+  std::vector<TreeDist::int16> a_in(a.n_splits);
+  std::vector<TreeDist::int16> b_in(b.n_splits);
+  for (TreeDist::int16 i = 0; i < a.n_splits; ++i) {
+    a_ptrs[i] = a.state[i];
+    a_in[i] = static_cast<TreeDist::int16>(a.in_split[i]);
+  }
+  for (TreeDist::int16 i = 0; i < b.n_splits; ++i) {
+    b_ptrs[i] = b.state[i];
+    b_in[i] = static_cast<TreeDist::int16>(b.in_split[i]);
+  }
+
+  return TreeDist::mutual_clustering_score(
+    a_ptrs.data(), a_in.data(), a.n_splits,
+    b_ptrs.data(), b_in.data(), b.n_splits,
+    a.n_bins, static_cast<TreeDist::int32>(n_tips),
+    scratch);
+}

tests/testthat/test-mci_impl.R

@@ -0,0 +1,53 @@
+test_that("Header mutual_clustering_score matches MutualClusteringInfo", {
+  skip_if_not_installed("TreeDist")
+  library(TreeTools)
+  library(TreeDist)
+
+  bal8 <- BalancedTree(8)
+  pec8 <- PectinateTree(8)
+  star8 <- StarTree(8)
+
+  tips <- TipLabels(bal8)
+  n_tip <- length(tips)
+  splits_bal <- as.Splits(bal8, tips)
+  splits_pec <- as.Splits(pec8, tips)
+  splits_star <- as.Splits(star8, tips)
+
+  # Score-only from the installable-header implementation
+  impl_score <- TreeDist:::cpp_mci_impl_score
+
+  impl_bal_pec <- impl_score(splits_bal, splits_pec, n_tip)
+  impl_bal_bal <- impl_score(splits_bal, splits_bal, n_tip)
+  impl_star <- impl_score(splits_bal, splits_star, n_tip)
+
+  # Reference from MutualClusteringInfo (unnormalized score)
+  ref_bal_pec <- MutualClusteringInfo(bal8, pec8)
+  ref_bal_bal <- MutualClusteringInfo(bal8, bal8)
+
+  expect_equal(impl_bal_pec, ref_bal_pec, tolerance = 1e-10)
+  expect_equal(impl_bal_bal, ref_bal_bal, tolerance = 1e-10)
+  expect_equal(impl_star, 0)
+})
+
+test_that("Header MCI covers exact-match early exit and partial LAP", {
+  skip_if_not_installed("TreeDist")
+  library(TreeTools)
+  library(TreeDist)
+  impl_score <- TreeDist:::cpp_mci_impl_score
+
+  # Two identical trees → all splits match exactly (early exit path)
+  bal20 <- BalancedTree(20)
+  tips <- TipLabels(bal20)
+  n_tip <- length(tips)
+  splits20 <- as.Splits(bal20, tips)
+
+  result <- impl_score(splits20, splits20, n_tip)
+  expect_equal(result, MutualClusteringInfo(bal20, bal20), tolerance = 1e-10)
+
+  # Trees that share some but not all splits → partial match + LAP
+  pec20 <- PectinateTree(20)
+  splits_pec20 <- as.Splits(pec20, tips)
+
+  result2 <- impl_score(splits20, splits_pec20, n_tip)
+  expect_equal(result2, MutualClusteringInfo(bal20, pec20), tolerance = 1e-10)
+})