diff --git a/src/ts_prune_reinsert.cpp b/src/ts_prune_reinsert.cpp index b2f99b8de..68a3ef092 100644 --- a/src/ts_prune_reinsert.cpp +++ b/src/ts_prune_reinsert.cpp @@ -323,11 +323,16 @@ DataSet build_reduced_dataset(const DataSet& ds, // Subset tip_states int tw = ds.total_words; red.tip_states.resize(static_cast(m) * tw); - for (int i = 0; i < m; ++i) { - int orig = tip_map[i]; - std::memcpy(&red.tip_states[static_cast(i) * tw], - &ds.tip_states[static_cast(orig) * tw], - tw * sizeof(uint64_t)); + // tw == 0 (zero Fitch blocks, e.g. every character constant/autapomorphic) + // leaves ds.tip_states empty; skip the memcpy rather than take the address + // of element 0 of an empty vector (undefined behaviour). + if (tw > 0) { + for (int i = 0; i < m; ++i) { + int orig = tip_map[i]; + std::memcpy(&red.tip_states[static_cast(i) * tw], + &ds.tip_states[static_cast(orig) * tw], + tw * sizeof(uint64_t)); + } } return red; @@ -428,18 +433,27 @@ void expand_and_reinsert( thread_local static int sa_delta_max = 0; #endif + // No Fitch words (e.g. every character constant/autapomorphic under equal + // weights): every insertion cost is identically zero. Skip the edge-set + // precompute and indirect length evaluation, which would otherwise take + // the address of element 0 of an empty vector (undefined behaviour; + // aborts under _GLIBCXX_ASSERTIONS) -- mirrors wagner_tree's guard. + const bool have_words = tw > 0; for (int tip : reinsert_order) { int new_internal = next_internal++; - const uint64_t* tip_prelim = - &ds.tip_states[static_cast(tip) * tw]; + const uint64_t* tip_prelim = have_words + ? &ds.tip_states[static_cast(tip) * tw] + : nullptr; // Exact insertion cost via directional edge sets: edge_set[D] = // combine(prelim[D], up[D]). Replaces the union-of-finals approximation // (final_[node] | final_[child]) that undercut insertion cost (~+30% Wagner // trees); mirrors the main Wagner builder. prelim is current here // (wagner_incremental_rescore maintains both prelim and final_). - compute_insertion_edge_sets(tree, ds, pr_edge_set, pr_up, pr_pre); + if (have_words) { + compute_insertion_edge_sets(tree, ds, pr_edge_set, pr_up, pr_pre); + } // Find best insertion edge via DFS from root int best_above = -1, best_below = -1; @@ -462,8 +476,11 @@ void expand_and_reinsert( if (lc < 0 || rc < 0) continue; // Evaluate edge (node, lc) - int extra = fitch_indirect_length_cached( - tip_prelim, &pr_edge_set[static_cast(lc) * tw], ds, best_extra); + int extra = have_words + ? fitch_indirect_length_cached( + tip_prelim, &pr_edge_set[static_cast(lc) * tw], ds, + best_extra) + : 0; if (extra < best_extra) { best_extra = extra; best_above = node; @@ -471,8 +488,11 @@ void expand_and_reinsert( } // Evaluate edge (node, rc) - extra = fitch_indirect_length_cached( - tip_prelim, &pr_edge_set[static_cast(rc) * tw], ds, best_extra); + extra = have_words + ? fitch_indirect_length_cached( + tip_prelim, &pr_edge_set[static_cast(rc) * tw], ds, + best_extra) + : 0; if (extra < best_extra) { best_extra = extra; best_above = node; @@ -496,7 +516,7 @@ void expand_and_reinsert( // is fully binary from root, so compute_insertion_edge_sets is exact and // safe. Tally Δ = exact_cost(E_bounded) − min_E exact_cost(E), per the // advisor: exact-suboptimality of the bounded choice, not raw edge flips. - if (best_below >= 0 && best_below != n_tip) { + if (have_words && best_below >= 0 && best_below != n_tip) { compute_insertion_edge_sets(tree, ds, sa_edge_set, sa_up, sa_pre); int exact_chosen = fitch_indirect_length_cached( tip_prelim, &sa_edge_set[static_cast(best_below) * tw], ds, INT_MAX); diff --git a/src/ts_sector.cpp b/src/ts_sector.cpp index c71b682b1..da2ac7466 100644 --- a/src/ts_sector.cpp +++ b/src/ts_sector.cpp @@ -72,8 +72,15 @@ static void compute_from_above_for_sector( // from_above[next] = fitch_join(from_above[node], prelim[sib]) // fitch_join: per-block, compute intersection; where empty, use union. - const uint64_t* sib_prelim = - &tree.prelim[static_cast(sib) * tw]; + // When tw == 0 (e.g. every character constant/autapomorphic under equal + // weights, leaving zero Fitch blocks), tree.prelim is empty and the + // block loop below never executes (ds.n_blocks == 0 too), so sib_prelim + // is never dereferenced -- but taking its address still needs guarding + // (address of element 0 of an empty vector is undefined behaviour; + // aborts under _GLIBCXX_ASSERTIONS). + const uint64_t* sib_prelim = (tw > 0) + ? &tree.prelim[static_cast(sib) * tw] + : nullptr; for (int b = 0; b < ds.n_blocks; ++b) { int off = ds.block_word_offset[b]; @@ -92,8 +99,10 @@ static void compute_from_above_for_sector( std::swap(from_above_cur, new_from_above); } - std::memcpy(from_above_out.data(), from_above_cur.data(), - tw * sizeof(uint64_t)); + if (tw > 0) { + std::memcpy(from_above_out.data(), from_above_cur.data(), + tw * sizeof(uint64_t)); + } } // ---- Conflict-guided sector selection ---- @@ -805,7 +814,11 @@ static ReducedDataset build_reduced_dataset_collapsed(const TreeState& tree, for (int i = 0; i < n_front; ++i) { // composite terminal states const int node = frontier[i]; const size_t dst = static_cast(i) * tw; - const uint64_t* src = (node < tree.n_tip) + // tw == 0 (zero Fitch blocks) leaves ds.tip_states/tree.prelim empty; + // guard the address-of (UB on an empty vector) -- the copy loop below + // is already a no-op (w < tw == 0), so src is never dereferenced. + const uint64_t* src = (tw == 0) ? nullptr + : (node < tree.n_tip) ? &ds.tip_states[static_cast(node) * tw] // real tip : &tree.prelim[static_cast(node) * tw]; // collapsed sub-clade for (int w = 0; w < tw; ++w) rd.data.tip_states[dst + w] = src[w]; @@ -921,12 +934,25 @@ static void build_ras_sector(ReducedDataset& rd, std::mt19937& rng) { // up-message buffer and preorder list each step. std::vector edge_set_up; std::vector edge_set_pre; + // When there are no Fitch words -- e.g. every character is constant or an + // autapomorphy, leaving the equal-weights dataset with zero blocks (the + // same all-uninformative case wagner_tree guards) -- rd.data.tip_states + // and edge_set are empty and every insertion cost is identically zero. + // Skip the edge-set precompute and the indirect length evaluation, which + // would otherwise take the address of element 0 of an empty vector + // (undefined behaviour; aborts under _GLIBCXX_ASSERTIONS/ASan). The DFS + // below still runs, so constraints are honoured and -- all costs being + // equal -- the first legal edge is chosen (mirrors wagner_tree's fix). + const bool have_words = rd.data.total_words > 0; for (int i = 2; i < n_real; ++i) { const int tip = order[i]; - const uint64_t* tip_prelim = - &rd.data.tip_states[static_cast(tip) * tw]; + const uint64_t* tip_prelim = have_words + ? &rd.data.tip_states[static_cast(tip) * tw] + : nullptr; - compute_insertion_edge_sets(t, rd.data, edge_set, edge_set_up, edge_set_pre); + if (have_words) { + compute_insertion_edge_sets(t, rd.data, edge_set, edge_set_up, edge_set_pre); + } int best_above = -1, best_below = -1, best_extra = INT_MAX; stack.clear(); @@ -939,14 +965,20 @@ static void build_ras_sector(ReducedDataset& rd, std::mt19937& rng) { int lc = t.left[ni]; int rc = t.right[ni]; if (lc >= 0) { - int extra = fitch_indirect_length_cached( - tip_prelim, &edge_set[static_cast(lc) * tw], rd.data, best_extra); + int extra = have_words + ? fitch_indirect_length_cached( + tip_prelim, &edge_set[static_cast(lc) * tw], + rd.data, best_extra) + : 0; if (extra < best_extra) { best_extra = extra; best_above = node; best_below = lc; } if (lc >= n_tip) stack.push_back(lc); } if (rc >= 0) { - int extra = fitch_indirect_length_cached( - tip_prelim, &edge_set[static_cast(rc) * tw], rd.data, best_extra); + int extra = have_words + ? fitch_indirect_length_cached( + tip_prelim, &edge_set[static_cast(rc) * tw], + rd.data, best_extra) + : 0; if (extra < best_extra) { best_extra = extra; best_above = node; best_below = rc; } if (rc >= n_tip) stack.push_back(rc); } diff --git a/src/ts_temper.cpp b/src/ts_temper.cpp index 41a903986..b995ae3d1 100644 --- a/src/ts_temper.cpp +++ b/src/ts_temper.cpp @@ -139,6 +139,11 @@ TemperResult stochastic_tbr_phase( double score = temper_full_rescore(tree, ds); double best_score = score; + // No informative characters: all trees have the same score. Skip the + // per-word clip/regraft scoring below, which would otherwise take the + // address of element 0 of the (empty, since total_words == 0) tree.prelim + // vector -- undefined behaviour (aborts under _GLIBCXX_ASSERTIONS). + if (ds.total_words == 0) return {best_score, score, 0, 0, 0}; const bool use_iw = std::isfinite(ds.concavity); const double eps = use_iw ? 1e-10 : 0.0; const double temperature = params.temperature; diff --git a/src/ts_tree.cpp b/src/ts_tree.cpp index 80357e218..a6a14aae4 100644 --- a/src/ts_tree.cpp +++ b/src/ts_tree.cpp @@ -52,9 +52,10 @@ void TreeState::load_tip_states(const DataSet& ds) { // T-262: bulk memcpy replaces per-element loop. Tip states occupy the // first n_tip * total_words entries of prelim/final_ (contiguous). size_t tip_bytes = static_cast(n_tip) * total_words * sizeof(uint64_t); - // All characters can be uninformative (e.g. every state a singleton), - // leaving zero blocks and a zero-length tip_states vector; its .data() - // is then permitted to be null, which memcpy's nonnull attribute forbids. + // All characters can be uninformative (e.g. every state constant or a + // singleton), leaving zero blocks and a zero-length tip_states vector; + // its .data() is then permitted to be null, which memcpy's nonnull + // attribute forbids. if (tip_bytes > 0) { std::memcpy(prelim.data(), ds.tip_states.data(), tip_bytes); std::memcpy(final_.data(), ds.tip_states.data(), tip_bytes); @@ -209,6 +210,14 @@ void TreeState::restore_prealloc_undo() { } void TreeState::save_node_state(int node) { + // No Fitch words (e.g. every character constant/autapomorphic under equal + // weights): there is no per-word state to save, so this is legitimately a + // no-op. Guard it explicitly -- the prealloc fast path below sizes its + // flat buffers to capacity * total_words, so total_words == 0 leaves them + // empty, and unconditionally memcpy'ing to/from element 0 of an empty + // vector is undefined behaviour (aborts under _GLIBCXX_ASSERTIONS). + if (total_words == 0) return; + // Fast path: use pre-allocated flat buffers (no heap allocation) if (prealloc_undo) { auto& u = *prealloc_undo; diff --git a/src/ts_wagner.cpp b/src/ts_wagner.cpp index 04c35cf6f..ada287f8b 100644 --- a/src/ts_wagner.cpp +++ b/src/ts_wagner.cpp @@ -611,6 +611,10 @@ std::vector wagner_goloboff_scores(const DataSet& ds) { int n_tip = ds.n_tips; int tw = ds.total_words; std::vector scores(n_tip, 0.0); + // No Fitch words (e.g. every character constant/autapomorphic): every tip + // is equally uninformative, so all scores are legitimately 0; skip the + // per-tip loop rather than take the address of an empty ds.tip_states. + if (tw == 0) return scores; for (int t = 0; t < n_tip; ++t) { double score = 0.0; @@ -644,6 +648,10 @@ std::vector wagner_entropy_scores(const DataSet& ds) { int n_tip = ds.n_tips; int tw = ds.total_words; std::vector scores(n_tip, 0.0); + // See wagner_goloboff_scores: no Fitch words means every score is + // legitimately 0; skip the loop rather than take the address of an empty + // ds.tip_states. + if (tw == 0) return scores; for (int t = 0; t < n_tip; ++t) { double score = 0.0; diff --git a/tests/testthat/test-ts-sector.R b/tests/testthat/test-ts-sector.R index f587bfdc9..892a00dc7 100644 --- a/tests/testthat/test-ts-sector.R +++ b/tests/testthat/test-ts-sector.R @@ -247,3 +247,84 @@ test_that("sector_diag with NA characters returns consistent scores", { expect_true(diag$clade_size >= 2) expect_true(diag$n_sector_tips == diag$clade_size + 1L) }) + +# ===== All-uninformative EW data: zero Fitch words (regression) ============= +# When every equal-weights character is constant or an autapomorphy (no state +# shared by >1 taxon), simplify_patterns removes every character, leaving +# DataSet::total_words == 0 and empty per-word state vectors (tip_states, +# tree.prelim, edge_set). Several sectorial code paths took the address of +# element 0 of these empty vectors -- undefined behaviour that aborts under +# hardened libstdc++ assertions / ASan -- including build_ras_sector() (the +# RAS-restart start-tree builder, only reached when rasStarts > 1) and +# compute_from_above_for_sector()/build_reduced_dataset_collapsed() (reached +# via ordinary rss_search/xss_search). 40 tips clears the sectorMinSize*2 +# gate (default sectorMinSize = 6) so a sector is always extracted. + +make_uninformative_ew <- function(n, n_char = 5L, seed) { + set.seed(seed) + mkchar <- function() { + states <- rep(1L, n) + # 6 tips get unique singleton states (2..7): no state is shared by more + # than one taxon, so every character is parsimony-uninformative. + states[sample(seq_len(n), 6L)] <- seq(2L, 7L) + states + } + mat <- vapply(seq_len(n_char), function(i) mkchar(), integer(n)) + rownames(mat) <- paste0("t", seq_len(n)) + MatrixToPhyDat(mat) +} + +test_that("Sectorial search handles all-uninformative EW data (zero Fitch words)", { + ds <- make_uninformative_ew(40L, seed = 1) + tree <- PectinateTree(names(ds)) + + res <- MaximizeParsimony(ds, tree = tree, maxReplicates = 2L, + targetHits = 1L, verbosity = 0L) + expect_s3_class(res[[1]], "phylo") + expect_equal(length(res[[1]]$tip.label), 40L) +}) + +test_that("build_ras_sector handles zero Fitch words (rasStarts > 1)", { + ds <- make_uninformative_ew(40L, seed = 1) + tree <- PectinateTree(names(ds)) + ctrl <- SearchControl(rasStarts = 3L) + + res <- MaximizeParsimony(ds, tree = tree, maxReplicates = 2L, + targetHits = 1L, verbosity = 0L, control = ctrl) + expect_s3_class(res[[1]], "phylo") + expect_equal(length(res[[1]]$tip.label), 40L) +}) + +test_that("build_reduced_dataset_collapsed handles zero Fitch words", { + ds <- make_uninformative_ew(40L, seed = 1) + tree <- PectinateTree(names(ds)) + ctrl <- SearchControl(sectorCollapseTarget = 6L) + + res <- MaximizeParsimony(ds, tree = tree, maxReplicates = 2L, + targetHits = 1L, verbosity = 0L, control = ctrl) + expect_s3_class(res[[1]], "phylo") + expect_equal(length(res[[1]]$tip.label), 40L) +}) + +test_that("expand_and_reinsert (prune-reinsert) handles zero Fitch words", { + ds <- make_uninformative_ew(40L, seed = 1) + tree <- PectinateTree(names(ds)) + ctrl <- SearchControl(pruneReinsertCycles = 1L) + + res <- MaximizeParsimony(ds, tree = tree, maxReplicates = 2L, + targetHits = 1L, verbosity = 0L, control = ctrl) + expect_s3_class(res[[1]], "phylo") + expect_equal(length(res[[1]]$tip.label), 40L) +}) + +test_that("stochastic_tbr_phase (annealing) handles zero Fitch words", { + ds <- make_uninformative_ew(40L, seed = 1) + tree <- PectinateTree(names(ds)) + ctrl <- SearchControl(annealCycles = 1L, annealPhases = 2L, + annealTStart = 5, annealTEnd = 0) + + res <- MaximizeParsimony(ds, tree = tree, maxReplicates = 2L, + targetHits = 1L, verbosity = 0L, control = ctrl) + expect_s3_class(res[[1]], "phylo") + expect_equal(length(res[[1]]$tip.label), 40L) +}) diff --git a/tests/testthat/test-ts-wagner.R b/tests/testthat/test-ts-wagner.R index d047b24cc..783e22c09 100644 --- a/tests/testthat/test-ts-wagner.R +++ b/tests/testthat/test-ts-wagner.R @@ -515,3 +515,34 @@ test_that("addition_order length/range/duplicate errors cleanly (T-323)", { "addition_order" ) }) + +# ===== All-uninformative data: zero Fitch words (regression) ================ +# When every character is constant or an autapomorphy (no state shared by +# more than one taxon), simplify_patterns removes every character, leaving +# DataSet::total_words == 0 and empty per-word state vectors. +# wagner_goloboff_scores()/wagner_entropy_scores() took the address of +# element 0 of the (empty) ds.tip_states vector -- undefined behaviour that +# aborts under hardened libstdc++ assertions / ASan. + +test_that("wagner_goloboff_scores/wagner_entropy_scores handle zero Fitch words", { + n <- 10L + set.seed(1) + # 3 characters, each a distinct permutation of 1..n: no state is shared by + # more than one taxon in any character, so every character is + # parsimony-uninformative (multi-character, avoiding the single-character + # vapply/t() degenerate case in prep_pd()). + mat <- vapply(seq_len(3L), function(i) as.character(sample(seq_len(n))), + character(n)) + rownames(mat) <- paste0("t", seq_len(n)) + pd <- MatrixToPhyDat(mat) + d <- prep_pd(pd) + + result <- TreeSearch:::ts_wagner_bias_bench( + d$contrast, d$tip_data, d$weight, d$levels, + min_steps = integer(0), concavity = -1, + bias = 1L, temperature = 1, n_reps = 1L, run_tbr = FALSE + ) + + expect_true(all(result$goloboff_scores == 0)) + expect_true(all(result$entropy_scores == 0)) +})