diff --git a/.AGENTS/memory/Optimizations.md b/.AGENTS/memory/Optimizations.md index 8f747bb2..aeaa5664 100644 --- a/.AGENTS/memory/Optimizations.md +++ b/.AGENTS/memory/Optimizations.md @@ -978,3 +978,12 @@ in TreeDist — every sort is `std::sort` / `std::stable_sort`. of those end below the cutoff (candidate rejected); the row bound exits ~2%, a column-max bound would add ~1%. Further gains need a tighter MCI upper bound or LAP warm starts between candidates. + +## Rogue SPR: parent dual prices (PR #215, 2026-10-14) +Full write-up: `dev/rogue_mci_bound_findings.md`. Parent column duals bound +82–88% of wasted LAPs; remaining LAPs warm-start (2–4×). −51% at 80 tips, +−65% at 150. Size-profile / compatibility bounds ~100× too loose; progressive +row-filling skips ≤1.5% of fill. Next lever: build child cost matrices +incrementally from the parent (97.6% of candidates drop one leaf); the fill is +now 37–55% of runtime. Parent prices are computed lazily (first candidate +past the row bound). diff --git a/DESCRIPTION b/DESCRIPTION index 0ab0027f..3e8aa811 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: TreeDist Title: Calculate and Map Distances Between Phylogenetic Trees -Version: 2.15.0.9006 +Version: 2.15.0.9007 Authors@R: c(person("Martin R.", "Smith", email = "martin.smith@durham.ac.uk", role = c("aut", "cre", "cph", "prg"), diff --git a/NEWS.md b/NEWS.md index fff6a3a0..ed8657f4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# TreeDist 2.15.0.9006 (development) +# TreeDist 2.15.0.9007 (development) ## Bug fixes @@ -8,6 +8,11 @@ ## Performance +- `SPRDist(method = "Rogue")` bounds each candidate's clustering information + with the optimal dual prices of the parent pair, skipping most assignment + problems and warm-starting the rest: c. 50% faster at 80 leaves and 65% at + 150 leaves on unrelated trees. + - Improve performance when comparing many similar trees with `ClusteringInfoDistance()`, `MatchingSplitDistance()`, `InfoRobinsonFoulds()` and `JaccardRobinsonFoulds()`. diff --git a/codemeta.json b/codemeta.json index dd718021..344d7c03 100644 --- a/codemeta.json +++ b/codemeta.json @@ -8,7 +8,7 @@ "codeRepository": "https://github.com/ms609/TreeDist/", "issueTracker": "https://github.com/ms609/TreeDist/issues/", "license": "https://spdx.org/licenses/GPL-3.0", - "version": "2.15.0.9006", + "version": "2.15.0.9007", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", @@ -571,7 +571,7 @@ "email": "martin.smith@durham.ac.uk" } ], - "name": "TreeDist: Distances between Phylogenetic Trees. R package version 2.15.0.9006", + "name": "TreeDist: Distances between Phylogenetic Trees. R package version 2.15.0.9007", "identifier": "10.32614/CRAN.package.TreeDist", "@id": "https://doi.org/10.32614/CRAN.package.TreeDist", "sameAs": "https://doi.org/10.32614/CRAN.package.TreeDist", diff --git a/dev/rogue_bound/add_mid.R b/dev/rogue_bound/add_mid.R new file mode 100644 index 00000000..e5e99fe3 --- /dev/null +++ b/dev/rogue_bound/add_mid.R @@ -0,0 +1,10 @@ +suppressMessages(library(TreeTools, lib.loc = "/root/Rlib")) +src <- readLines("make_workload.R") +end <- grep("^set.seed", src)[1] - 1 +eval(parse(text = src[1:end])) +wl <- readRDS("workloads.rds") +set.seed(777) +wl$mid80 <- lapply(1:12, function(i) NearPair(80, 15)) +wl$mid150 <- lapply(1:6, function(i) NearPair(150, 30)) +saveRDS(wl, "workloads.rds") +cat(names(wl), "\n") diff --git a/dev/rogue_bound/alt_summary.R b/dev/rogue_bound/alt_summary.R new file mode 100644 index 00000000..04b79f76 --- /dev/null +++ b/dev/rogue_bound/alt_summary.R @@ -0,0 +1,17 @@ +dir <- "/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/csv2" +out <- NULL +for (f in c("rand30", "rand30d2", "rand80", "rand150", "near80", "near150")) { + x <- read.csv(file.path(dir, paste0(f, ".csv"))) + lapCalls <- x$ub_row >= x$min_score + w <- x[lapCalls & x$lap_exact < x$min_score, ] + Rej <- function(b) { v <- w[[b]]; v[is.na(v)] <- Inf; 100 * mean(v < w$min_score) } + out <- rbind(out, data.frame( + wl = f, n_wasted = nrow(w), + cr = Rej("ub_cr"), rc = Rej("ub_rc"), best_cr_rc = Rej("ub_alt1"), + alt3 = Rej("ub_alt3"), alt10 = Rej("ub_alt10"), + pv = Rej("ub_pv_minv"), pv_max = Rej("ub_pv_maxv"), pv_mean = Rej("ub_pv_mean"), + pu = Rej("ub_pu"), pv1 = Rej("ub_pv_alt1"), pv3 = Rej("ub_pv_alt3"), + pvarr = Rej("ub_pv_warmarr"), prev = Rej("ub_prev"), + arr = Rej("ub_arr"), rt = Rej("ub_rt"))) +} +print(format(out, digits = 3), row.names = FALSE) diff --git a/dev/rogue_bound/ana/ana b/dev/rogue_bound/ana/ana new file mode 100644 index 00000000..7591b059 Binary files /dev/null and b/dev/rogue_bound/ana/ana differ diff --git a/dev/rogue_bound/ana/ana.cpp b/dev/rogue_bound/ana/ana.cpp new file mode 100644 index 00000000..291a332b --- /dev/null +++ b/dev/rogue_bound/ana/ana.cpp @@ -0,0 +1,682 @@ +// Offline evaluation of upper bounds on the bounded MCI calls dumped by the +// instrumented TreeDist build (spr.cpp / instr_dump()). +// +// Usage: ana +#include +static unsigned long long g_init[2], g_arr[2], g_aug[2], g_n[2]; +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +#include "lap_hooked.h" + +using namespace TreeDist; +using i128 = __int128; + +static double lg2tab[4096]; +static void InitLg2() { + lg2tab[0] = 0; + for (int i = 1; i < 4096; ++i) lg2tab[i] = std::log2(double(i)); +} + +struct SplitSet { + int n_tip = 0, n_split = 0, nw = 0; + std::vector bits; + std::vector in; + const uint64_t *row(int i) const { return &bits[size_t(i) * nw]; } +}; + +static bool ReadSplits(FILE *f, int n_tip, SplitSet &s) { + int hdr[2]; + if (fread(hdr, sizeof(int), 2, f) != 2) return false; + s.n_tip = n_tip; + s.n_split = hdr[0]; + const int nc = hdr[1]; + s.nw = (n_tip + 63) / 64; + s.bits.assign(size_t(s.n_split) * s.nw, 0); + s.in.assign(s.n_split, 0); + std::vector buf(nc); + for (int r = 0; r != s.n_split; ++r) { + if (fread(buf.data(), 1, nc, f) != size_t(nc)) return false; + for (int c = 0; c != nc; ++c) { + s.bits[size_t(r) * s.nw + c / 8] |= uint64_t(buf[c]) << (8 * (c % 8)); + } + int cnt = 0; + for (int w = 0; w != s.nw; ++w) cnt += __builtin_popcountll(s.bits[size_t(r) * s.nw + w]); + s.in[r] = cnt; + } + return true; +} + +// ---- Weights exactly as mutual_clustering() ---- +static inline double IcMatching(int a, int b, int n) { + return (a + b) * lg2tab[n] - a * lg2tab[a] - b * lg2tab[b]; +} + +static inline cost CellCost(int na, int nb, int a_and_b, int n) { + const int nA = n - na, nB = n - nb; + const int a_and_B = na - a_and_b; + const int A_and_b = nb - a_and_b; + const int A_and_B = nA - A_and_b; + if (a_and_b == A_and_b && a_and_b == a_and_B && a_and_b == A_and_B) return BIG; + const double lg2_n = lg2tab[n]; + const double offset_a = lg2_n - lg2tab[na]; + const double offset_A = lg2_n - lg2tab[nA]; + const double lg2_nb = lg2tab[nb], lg2_nB = lg2tab[nB]; + const double ic_sum = + a_and_b * (lg2tab[a_and_b] + offset_a - lg2_nb) + + a_and_B * (lg2tab[a_and_B] + offset_a - lg2_nB) + + A_and_b * (lg2tab[A_and_b] + offset_A - lg2_nb) + + A_and_B * (lg2tab[A_and_B] + offset_A - lg2_nB); + const double max_over_tips = double(BIG) / n; + return BIG - cost(ic_sum * max_over_tips); +} + +// Reduced square LAP problem after removing exact matches. +struct Problem { + int n_tip = 0, d = 0; + double exact_part = 0; // exact-match score, already divided by n + std::vector rows, cols; // original split indices; -1 = padding + CostMatrix C; + double entropy_bound = 0; +}; + +static void Build(const SplitSet &a, const SplitSet &b, Problem &p) { + const int n = a.n_tip; + p.n_tip = n; + std::vector a_match(a.n_split, -1), b_match(b.n_split, -1); + double exact_ic = 0; + double ent_a = 0, ent_b = 0; + for (int i = 0; i < a.n_split; ++i) ent_a += IcMatching(a.in[i], n - a.in[i], n); + for (int j = 0; j < b.n_split; ++j) ent_b += IcMatching(b.in[j], n - b.in[j], n); + p.entropy_bound = std::min(ent_a, ent_b) / n; + for (int i = 0; i < a.n_split; ++i) { + const int na = a.in[i], nA = n - na; + for (int j = 0; j < b.n_split; ++j) { + int ab = 0; + for (int w = 0; w < a.nw; ++w) ab += __builtin_popcountll(a.row(i)[w] & b.row(j)[w]); + const int nb = b.in[j]; + const int aB = na - ab, Ab = nb - ab, AB = nA - Ab; + if ((!aB && !Ab) || (!ab && !AB)) { + exact_ic += IcMatching(na, nA, n); + a_match[i] = j; + b_match[j] = i; + break; + } + } + } + p.exact_part = exact_ic / n; + p.rows.clear(); + p.cols.clear(); + for (int i = 0; i < a.n_split; ++i) if (a_match[i] < 0) p.rows.push_back(i); + for (int j = 0; j < b.n_split; ++j) if (b_match[j] < 0) p.cols.push_back(j); + p.d = std::max(p.rows.size(), p.cols.size()); + while (int(p.rows.size()) < p.d) p.rows.push_back(-1); + while (int(p.cols.size()) < p.d) p.cols.push_back(-1); + p.C.resize(p.d); + for (int r = 0; r < p.d; ++r) { + const int i = p.rows[r]; + for (int c = 0; c < p.d; ++c) { + const int j = p.cols[c]; + if (i < 0 || j < 0) { p.C(r, c) = BIG; continue; } + int ab = 0; + for (int w = 0; w < a.nw; ++w) ab += __builtin_popcountll(a.row(i)[w] & b.row(j)[w]); + p.C(r, c) = CellCost(a.in[i], b.in[j], ab, n); + } + } +} + +static inline double ToScore(const Problem &p, i128 lb) { + // Same arithmetic as mutual_clustering(): (BIG * dim - cost) / BIG + return double(cost(i128(BIG) * p.d - lb)) * (1.0 / double(BIG)) + p.exact_part; +} + +// LB(v) = sum v + sum_i min_j (c_ij - v_j): valid lower bound on LAP cost. +static i128 LowerBound(const Problem &p, const std::vector &v, + std::vector *u_out = nullptr) { + i128 lb = 0; + for (int j = 0; j < p.d; ++j) lb += v[j]; + if (u_out) u_out->resize(p.d); + for (int i = 0; i < p.d; ++i) { + const cost *row = p.C.row(i); + cost m = std::numeric_limits::max(); + for (int j = 0; j < p.d; ++j) m = std::min(m, row[j] - v[j]); + lb += m; + if (u_out) (*u_out)[i] = m; + } + return lb; +} + +// From row duals u: v_j = min_i (c_ij - u_i). LB = sum u + sum v. +static i128 LowerBoundU(const Problem &p, const std::vector &u, + std::vector *v_out = nullptr) { + i128 lb = 0; + for (int i = 0; i < p.d; ++i) lb += u[i]; + std::vector v(p.d, std::numeric_limits::max()); + for (int i = 0; i < p.d; ++i) { + const cost *row = p.C.row(i); + for (int j = 0; j < p.d; ++j) v[j] = std::min(v[j], row[j] - u[i]); + } + for (int j = 0; j < p.d; ++j) lb += v[j]; + if (v_out) *v_out = v; + return lb; +} + +// Alternating (block coordinate) dual ascent from v, k rounds; returns LB. +static i128 Alternate(const Problem &p, std::vector v, int k) { + std::vector u; + i128 lb = LowerBound(p, v, &u); + for (int r = 0; r < k; ++r) { + lb = std::max(lb, LowerBoundU(p, u, &v)); + lb = std::max(lb, LowerBound(p, v, &u)); + } + return lb; +} + +struct LapResult { + cost opt = 0; + std::vector v; // final column duals + std::vector rowsol; + ana::LapWork work; + uint64_t ticks = 0; + // Augmentation-phase bound: augmentations done when LB first rejects + int aug_total = 0, aug_needed = -1; + double ub_rt = 0, ub_arr = 0, ub_col = 0; + int ub_row_dummy = 0; +}; + +static LapResult Solve(Problem &p, const cost *v_init, bool track, + double min_score) { + LapResult res; + std::vector rowsol(p.d); + std::vector colsol(p.d); + LapScratch scratch; + CostMatrix work_c = p.C; // lap may transpose + if (!track) { + const uint64_t t0 = __rdtsc(); + res.opt = ana::lap_hooked(p.d, work_c, rowsol, colsol, scratch, v_init, + [](int, const std::vector &, int, int) {}, + res.work); + res.ticks = __rdtsc() - t0; + } else { + std::vector vcopy; + auto hook = [&](int phase, const std::vector &v, int done, int total) { + vcopy.assign(v.begin(), v.begin() + p.d); + if (phase == 0) res.ub_col = ToScore(p, LowerBound(p, vcopy)); + else if (phase == 1) res.ub_rt = ToScore(p, LowerBound(p, vcopy)); + else if (phase == 2) { + res.ub_arr = ToScore(p, LowerBound(p, vcopy)); + res.aug_total = total; + if (res.ub_arr < min_score) res.aug_needed = 0; + } else if (phase == 3 && res.aug_needed < 0) { + // Check at ~10% steps to limit analysis cost + const int step = std::max(1, total / 10); + if (done % step == 0 || done == total) { + if (ToScore(p, LowerBound(p, vcopy)) < min_score) res.aug_needed = done; + } + } + }; + res.opt = ana::lap_hooked(p.d, work_c, rowsol, colsol, scratch, v_init, + hook, res.work); + } + res.v.assign(scratch.v.begin(), scratch.v.begin() + p.d); + res.rowsol.assign(rowsol.begin(), rowsol.end()); + return res; +} + +// ---- Mapping splits between leaf sets ---- +using Key = std::vector; +struct KeyHash { + size_t operator()(const Key &k) const { + uint64_t h = 1469598103934665603ULL; + for (auto w : k) { h ^= w; h *= 1099511628211ULL; h ^= h >> 29; } + return h; + } +}; + +// A set of splits on parent leaf indexing with associated dual values. +struct PriceSource { + int nwp = 0; // words in parent indexing + std::vector sets; + std::vector val; +}; + +static Key ToParent(const SplitSet &s, int i, const std::vector &map, int nwp) { + Key k(nwp, 0); + for (int t = 0; t < s.n_tip; ++t) { + if ((s.row(i)[t / 64] >> (t % 64)) & 1) { + const int pt = map[t]; + k[pt / 64] |= uint64_t(1) << (pt % 64); + } + } + return k; +} + +// Restrict to leaf set L, polarized so that reference leaf r is outside. +static Key Restrict(const Key &s, const Key &L, int r) { + Key k(s.size()); + const bool flip = (s[r / 64] >> (r % 64)) & 1; + for (size_t w = 0; w < s.size(); ++w) k[w] = (flip ? ~s[w] : s[w]) & L[w]; + return k; +} + +enum Agg { AGG_MIN_V, AGG_MAX_V, AGG_MEAN_V }; + +// Map source prices onto the target problem's columns (if cols) or rows. +// Returns per-index value; `mapped` flags which received one. +static void MapPrices(const PriceSource &src, const SplitSet &target, + const std::vector &target_idx, + const std::vector &map, Agg agg, + std::vector &out, std::vector &mapped, + int &n_multi) { + const int nwp = src.nwp; + Key L(nwp, 0); + for (int t = 0; t < target.n_tip; ++t) L[map[t] / 64] |= uint64_t(1) << (map[t] % 64); + const int r = map[0]; + std::unordered_map, KeyHash> table; + for (size_t q = 0; q < src.sets.size(); ++q) { + table[Restrict(src.sets[q], L, r)].push_back(int(q)); + } + const int d = target_idx.size(); + out.assign(d, 0); + mapped.assign(d, 0); + n_multi = 0; + for (int c = 0; c < d; ++c) { + const int j = target_idx[c]; + if (j < 0) continue; + const Key k = Restrict(ToParent(target, j, map, nwp), L, r); + auto it = table.find(k); + if (it == table.end()) continue; + const auto &qs = it->second; + if (qs.size() > 1) ++n_multi; + cost acc; + if (agg == AGG_MIN_V) { + acc = std::numeric_limits::max(); + for (int q : qs) acc = std::min(acc, src.val[q]); + } else if (agg == AGG_MAX_V) { + acc = std::numeric_limits::min(); + for (int q : qs) acc = std::max(acc, src.val[q]); + } else { + long double s = 0; + for (int q : qs) s += src.val[q]; + acc = cost(s / qs.size()); + } + out[c] = acc; + mapped[c] = 1; + } +} + +// Fill unmapped column prices with the best feasible value given rows' +// minima over mapped columns: v_j = min_i (c_ij - u_i). +static bool CompleteV(const Problem &p, std::vector &v, + const std::vector &mapped) { + int nm = 0; + for (char m : mapped) nm += m; + if (nm == 0) return false; + if (nm == p.d) return true; + std::vector u(p.d); + for (int i = 0; i < p.d; ++i) { + cost m = std::numeric_limits::max(); + for (int j = 0; j < p.d; ++j) if (mapped[j]) m = std::min(m, p.C(i, j) - v[j]); + u[i] = m; + } + for (int j = 0; j < p.d; ++j) { + if (mapped[j]) continue; + cost m = std::numeric_limits::max(); + for (int i = 0; i < p.d; ++i) m = std::min(m, p.C(i, j) - u[i]); + v[j] = m; + } + return true; +} + +static bool CompleteU(const Problem &p, std::vector &u, + const std::vector &mapped) { + int nm = 0; + for (char m : mapped) nm += m; + if (nm == 0) return false; + if (nm == p.d) return true; + std::vector v(p.d); + for (int j = 0; j < p.d; ++j) { + cost m = std::numeric_limits::max(); + for (int i = 0; i < p.d; ++i) if (mapped[i]) m = std::min(m, p.C(i, j) - u[i]); + v[j] = m; + } + for (int i = 0; i < p.d; ++i) { + if (mapped[i]) continue; + cost m = std::numeric_limits::max(); + for (int j = 0; j < p.d; ++j) m = std::min(m, p.C(i, j) - v[j]); + u[i] = m; + } + return true; +} + +// Size-only bound: max-weight matching where each weight is the largest +// MCI attainable by two splits of these sizes (nested configurations). +static double SizeBound(const Problem &p, const SplitSet &a, const SplitSet &b) { + Problem q; + q.n_tip = p.n_tip; q.d = p.d; q.exact_part = p.exact_part; + q.C.resize(p.d); + const int n = p.n_tip; + for (int r = 0; r < p.d; ++r) { + for (int c = 0; c < p.d; ++c) { + const int i = p.rows[r], j = p.cols[c]; + if (i < 0 || j < 0) { q.C(r, c) = BIG; continue; } + const int na = a.in[i], nb = b.in[j]; + const cost c1 = CellCost(na, nb, std::min(na, nb), n); + const cost c2 = CellCost(na, nb, std::max(0, na + nb - n), n); + q.C(r, c) = std::min(c1, c2); + } + } + LapResult res = Solve(q, nullptr, false, 0); + return ToScore(q, res.opt); +} + +// Progressive priced bound: rows filled in `order`; unfilled rows use the +// provable bound c_ij >= BIG - floor(BIG * min(h_i, h_j)) - 1. Returns the +// fraction of rows that must be filled before the bound rejects (2 = never). +static double ProgressiveFraction(const Problem &p, const SplitSet &a, + const SplitSet &b, const std::vector &v, + double min_score, int order_kind) { + const int n = p.n_tip, d = p.d; + std::vector hr(d, 0), hc(d, 0); + for (int r = 0; r < d; ++r) if (p.rows[r] >= 0) hr[r] = IcMatching(a.in[p.rows[r]], n - a.in[p.rows[r]], n) / n; + for (int c = 0; c < d; ++c) if (p.cols[c] >= 0) hc[c] = IcMatching(b.in[p.cols[c]], n - b.in[p.cols[c]], n) / n; + std::vector exact_u(d), loose_u(d); + i128 lb = 0; + for (int j = 0; j < d; ++j) lb += v[j]; + for (int i = 0; i < d; ++i) { + cost mu = std::numeric_limits::max(), ml = std::numeric_limits::max(); + for (int j = 0; j < d; ++j) { + mu = std::min(mu, p.C(i, j) - v[j]); + const cost lower = BIG - cost(std::floor(double(BIG) * std::min(hr[i], hc[j]))) - 1; + ml = std::min(ml, std::max(cost(0), lower) - v[j]); + } + exact_u[i] = mu; + loose_u[i] = ml; + lb += ml; + } + std::vector ord(d); + for (int i = 0; i < d; ++i) ord[i] = i; + if (order_kind == 1) { + // Oracle-free heuristic: largest provable slack first + std::sort(ord.begin(), ord.end(), [&](int x, int y) { + return hr[x] > hr[y]; + }); + } else if (order_kind == 2) { + // Oracle: largest actual tightening first (best possible order) + std::sort(ord.begin(), ord.end(), [&](int x, int y) { + return exact_u[x] - loose_u[x] > exact_u[y] - loose_u[y]; + }); + } + if (ToScore(p, lb) < min_score) return 0; + for (int t = 0; t < d; ++t) { + const int i = ord[t]; + lb += exact_u[i] - loose_u[i]; + if (ToScore(p, lb) < min_score) return double(t + 1) / d; + } + return 2; +} + +struct Child { + int cand, kind, n_tip; + std::vector map; + double min_score, exact, bounded, info; + SplitSet a, b; +}; + +int main(int argc, char **argv) { + if (argc < 3) { fprintf(stderr, "usage\n"); return 1; } + InitLg2(); + FILE *f = fopen(argv[1], "rb"); + FILE *out = fopen(argv[2], "w"); + const bool do_size = argc > 3 && std::string(argv[3]) == "size"; + fprintf(out, "parent,cand,kind,n_tip,d,exact,lap_exact,min_score," + "ub_entropy,ub_row,ub_cr,ub_rc,ub_rt,ub_arr," + "ub_alt1,ub_alt3,ub_alt10," + "ub_pv_minv,ub_pv_maxv,ub_pv_mean,ub_pu,ub_pv_alt1,ub_pv_alt3," + "ub_pv_warmarr,ub_prev,ub_size,frac_mapped,n_multi," + "aug_total,aug_needed," + "cold_free_init,cold_free_arr,cold_arr,cold_aug,cold_ticks," + "warm_free_init,warm_free_arr,warm_arr,warm_aug,warm_ticks," + "warm_ok,prev_ticks,parent_n_tip,prog_nat,prog_h,prog_oracle\n"); + int parent_id = -1; + int parent_n_tip = 0; + PriceSource parent_v, parent_u; + PriceSource prev_v; + bool have_prev = false; + int prev_parent = -2; + long n_records = 0, n_mismatch = 0; + int tag; + while (fread(&tag, sizeof(int), 1, f) == 1) { + if (tag == 1) { + int n_tip; + double info; + fread(&n_tip, sizeof(int), 1, f); + fread(&info, sizeof(double), 1, f); + SplitSet A, B; + ReadSplits(f, n_tip, A); + ReadSplits(f, n_tip, B); + ++parent_id; + parent_n_tip = n_tip; + have_prev = false; + Problem p; + Build(A, B, p); + const int nwp = (n_tip + 63) / 64; + std::vector ident(n_tip); + for (int t = 0; t < n_tip; ++t) ident[t] = t; + parent_v = PriceSource{nwp, {}, {}}; + parent_u = PriceSource{nwp, {}, {}}; + if (p.d > 1) { + LapResult res = Solve(p, nullptr, false, 0); + // Check optimality of the final duals + std::vector u; + const i128 lb = LowerBound(p, res.v, &u); + if (lb != i128(res.opt)) { + fprintf(stderr, "parent %d: dual gap %g\n", parent_id, double(i128(res.opt) - lb)); + } + { + cost mn = *std::min_element(res.v.begin(), res.v.end()); + cost mx = *std::max_element(res.v.begin(), res.v.end()); + static double worst = 0; + const double spread = double(mx - mn) / double(BIG); + if (spread > worst) { worst = spread; fprintf(stderr, "parent %d d=%d: dual spread %.3f BIG\n", parent_id, p.d, spread); } + } + for (int c = 0; c < p.d; ++c) { + if (p.cols[c] < 0) continue; + parent_v.sets.push_back(ToParent(B, p.cols[c], ident, nwp)); + parent_v.val.push_back(res.v[c]); + } + for (int r = 0; r < p.d; ++r) { + if (p.rows[r] < 0) continue; + parent_u.sets.push_back(ToParent(A, p.rows[r], ident, nwp)); + parent_u.val.push_back(u[r]); + } + } + continue; + } + if (tag != 2) { fprintf(stderr, "bad tag %d\n", tag); return 1; } + Child ch; + int hdr[3]; + fread(hdr, sizeof(int), 3, f); + ch.cand = hdr[0]; ch.kind = hdr[1]; ch.n_tip = hdr[2]; + ch.map.resize(ch.n_tip); + fread(ch.map.data(), sizeof(int), ch.n_tip, f); + double vals[4]; + fread(vals, sizeof(double), 4, f); + ch.min_score = vals[0]; ch.exact = vals[1]; ch.bounded = vals[2]; ch.info = vals[3]; + ReadSplits(f, ch.n_tip, ch.a); + ReadSplits(f, ch.n_tip, ch.b); + ++n_records; + + Problem p; + Build(ch.a, ch.b, p); + if (p.d == 0 || ch.a.n_split == 0 || ch.b.n_split == 0 || + int(p.rows.size()) == 0) { + continue; + } + // All-exact case + bool all_exact = true; + for (int r = 0; r < p.d; ++r) if (p.rows[r] >= 0 && p.cols[r] >= 0) all_exact = false; + if (p.d == 0 || all_exact) continue; + + LapResult cold = Solve(p, nullptr, true, ch.min_score); + const double lap_exact = ToScore(p, cold.opt); + if (std::fabs(lap_exact - ch.exact) > 1e-9) ++n_mismatch; + + // Row and simple bounds + std::vector zero(p.d, 0), u0; + const double ub_row = ToScore(p, LowerBound(p, zero, &u0)); + std::vector v_rc; + const double ub_rc = ToScore(p, LowerBoundU(p, u0, &v_rc)); + // column-first then alternating + std::vector vcol(p.d, std::numeric_limits::max()); + for (int i = 0; i < p.d; ++i) for (int j = 0; j < p.d; ++j) vcol[j] = std::min(vcol[j], p.C(i, j)); + const double ub_alt1 = ToScore(p, std::max(Alternate(p, vcol, 1), LowerBoundU(p, u0))); + const double ub_alt3 = ToScore(p, Alternate(p, vcol, 3)); + const double ub_alt10 = ToScore(p, Alternate(p, vcol, 10)); + + // Parent prices + double ub_pv[3] = {NAN, NAN, NAN}, ub_pu = NAN, ub_pv_alt1 = NAN, ub_pv_alt3 = NAN; + double ub_warmarr = NAN, frac_mapped = NAN; + int n_multi = 0; + LapResult warm; + bool warm_ok = false; + std::vector v_best; + if (!parent_v.sets.empty()) { + for (int a = 0; a < 3; ++a) { + std::vector v; + std::vector mapped; + int nm; + MapPrices(parent_v, ch.b, p.cols, ch.map, Agg(a), v, mapped, nm); + if (a == 0) { + int cnt = 0; + for (char m : mapped) cnt += m; + frac_mapped = double(cnt) / p.d; + n_multi = nm; + } + if (CompleteV(p, v, mapped)) { + ub_pv[a] = ToScore(p, LowerBound(p, v)); + if (a == 0) v_best = v; + } + } + { + std::vector u; + std::vector mapped; + int nm; + MapPrices(parent_u, ch.a, p.rows, ch.map, AGG_MAX_V, u, mapped, nm); + if (CompleteU(p, u, mapped)) ub_pu = ToScore(p, LowerBoundU(p, u)); + } + if (!v_best.empty()) { + ub_pv_alt1 = ToScore(p, Alternate(p, v_best, 1)); + ub_pv_alt3 = ToScore(p, Alternate(p, v_best, 3)); + // Warm-started JV from parent prices; bound after ARR + double ub_after_arr = NAN; + ana::LapWork w; + LapResult wr; + { + std::vector rowsol(p.d); + std::vector colsol(p.d); + LapScratch scratch; + CostMatrix work_c = p.C; + std::vector vv; + auto hook = [&](int phase, const std::vector &v, int, int) { + if (phase == 2) { + vv.assign(v.begin(), v.begin() + p.d); + ub_after_arr = ToScore(p, LowerBound(p, vv)); + } + }; + wr.opt = ana::lap_hooked(p.d, work_c, rowsol, colsol, scratch, + v_best.data(), hook, wr.work); + } + ub_warmarr = ub_after_arr; + // Timed warm run without hooks + warm = Solve(p, v_best.data(), false, 0); + warm_ok = warm.opt == cold.opt; + if (!warm_ok) { + const i128 lbw = LowerBound(p, warm.v), lbc = LowerBound(p, cold.v); + fprintf(stderr, "mismatch d=%d warm-cold %lld; warm gap %lld; cold gap %lld (units; 1 bit = %lld)\n", + p.d, (long long)(warm.opt - cold.opt), (long long)(warm.opt - lbw), + (long long)(cold.opt - lbc), (long long)BIG); + } + } + } + // Previous call's optimal prices (same parent) + double ub_prev = NAN; + unsigned long long prev_ticks = 0; + if (have_prev && prev_parent == parent_id) { + std::vector v; + std::vector mapped; + int nm; + MapPrices(prev_v, ch.b, p.cols, ch.map, AGG_MIN_V, v, mapped, nm); + if (CompleteV(p, v, mapped)) { + ub_prev = ToScore(p, LowerBound(p, v)); + LapResult pr = Solve(p, v.data(), false, 0); + prev_ticks = pr.ticks; + } + } + // Remember this call's optimal duals + { + const int nwp = parent_v.nwp ? parent_v.nwp : (ch.n_tip + 63) / 64; + prev_v = PriceSource{nwp, {}, {}}; + for (int c = 0; c < p.d; ++c) { + if (p.cols[c] < 0) continue; + prev_v.sets.push_back(ToParent(ch.b, p.cols[c], ch.map, nwp)); + prev_v.val.push_back(cold.v[c]); + } + have_prev = true; + prev_parent = parent_id; + } + const double ub_size = do_size ? SizeBound(p, ch.a, ch.b) : NAN; + double prog[3] = {NAN, NAN, NAN}; + if (!v_best.empty()) { + for (int k = 0; k < 3; ++k) prog[k] = ProgressiveFraction(p, ch.a, ch.b, v_best, ch.min_score, k); + } + // Timed cold run without hooks + LapResult cold_t = Solve(p, nullptr, false, 0); + { + const int kk = (cold.ub_row_dummy, 0); + (void)kk; + g_init[0] += cold_t.work.t_init; g_arr[0] += cold_t.work.t_arr; g_aug[0] += cold_t.work.t_aug; ++g_n[0]; + g_init[1] += warm.work.t_init; g_arr[1] += warm.work.t_arr; g_aug[1] += warm.work.t_aug; ++g_n[1]; + } + + fprintf(out, "%d,%d,%d,%d,%d,%.12g,%.12g,%.12g," + "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g," + "%.12g,%.12g,%.12g," + "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g," + "%.12g,%.12g,%.12g,%.4g,%d," + "%d,%d," + "%lld,%lld,%lld,%lld,%llu," + "%lld,%lld,%lld,%lld,%llu,%d,%llu,%d,%.4g,%.4g,%.4g\n", + parent_id, ch.cand, ch.kind, ch.n_tip, p.d, ch.exact, lap_exact, + ch.min_score, p.entropy_bound, ub_row, cold.ub_col, ub_rc, + cold.ub_rt, cold.ub_arr, ub_alt1, ub_alt3, ub_alt10, + ub_pv[0], ub_pv[1], ub_pv[2], ub_pu, ub_pv_alt1, ub_pv_alt3, + ub_warmarr, ub_prev, ub_size, frac_mapped, n_multi, + cold.aug_total, cold.aug_needed, + cold_t.work.n_free_after_init, cold_t.work.n_free_after_arr, + cold_t.work.arr_steps, cold_t.work.aug_inner + cold_t.work.aug_scan, + (unsigned long long)cold_t.ticks, + warm.work.n_free_after_init, warm.work.n_free_after_arr, + warm.work.arr_steps, warm.work.aug_inner + warm.work.aug_scan, + (unsigned long long)warm.ticks, int(warm_ok), prev_ticks, parent_n_tip, + prog[0], prog[1], prog[2]); + } + fprintf(stderr, "records %ld, exact mismatches %ld\n", n_records, n_mismatch); + for (int k = 0; k < 2; ++k) { + const double tot = double(g_init[k] + g_arr[k] + g_aug[k]); + fprintf(stderr, "%s LAP: mean %.1f us; init %.1f%%, ARR %.1f%%, augment %.1f%%\n", k ? "warm" : "cold", + tot / g_n[k] / 2100.0, 100 * g_init[k] / tot, 100 * g_arr[k] / tot, 100 * g_aug[k] / tot); + } + fclose(out); + fclose(f); +} diff --git a/dev/rogue_bound/ana/cost_matrix.h b/dev/rogue_bound/ana/cost_matrix.h new file mode 100644 index 00000000..33a48f91 --- /dev/null +++ b/dev/rogue_bound/ana/cost_matrix.h @@ -0,0 +1,277 @@ +#ifndef TREEDIST_COST_MATRIX_H_ +#define TREEDIST_COST_MATRIX_H_ + +// TreeDist CostMatrix — cache-aligned, block-padded cost matrix for the LAP. +// +// Rcpp-free public API. The Rcpp-dependent constructor +// (CostMatrix(Rcpp::NumericMatrix)) lives in TreeDist's own src/lap.h. + +#include "types.h" + +#include +#include +#include +#include + +namespace TreeDist { + +class CostMatrix { +private: + std::size_t dim_; + std::size_t dim8_; + alignas(64) std::vector data_; + alignas(64) std::vector t_data_; + bool transposed_; + + static std::size_t block_containing(const std::size_t x) { + return ((x + BLOCK_SIZE - 1) / BLOCK_SIZE) * BLOCK_SIZE; + } + +public: + // Default constructor for pooled instances (zero-size, no allocation). + CostMatrix() : dim_(0), dim8_(0), transposed_(false) {} + + explicit CostMatrix(std::size_t dim) + : dim_(dim), + dim8_(block_containing(dim_)), + data_(std::vector(dim8_ * dim8_)), + t_data_(std::vector(dim8_ * dim8_)), + transposed_(false) {} + + // Resize for reuse. Only reallocates when the new dimension exceeds the + // current buffer capacity; otherwise just updates dim_/dim8_ and marks + // the transpose as stale. + void resize(std::size_t new_dim) { + dim_ = new_dim; + dim8_ = block_containing(new_dim); + const std::size_t needed = dim8_ * dim8_; + if (data_.size() < needed) { + data_.resize(needed); + t_data_.resize(needed); + } + transposed_ = false; + } + + void reset() noexcept { transposed_ = false; } + + [[nodiscard]] std::size_t dim() const noexcept { return dim_; } + [[nodiscard]] std::size_t dim8() const noexcept { return dim8_; } + + // ---- Element access ---- + + cost& operator()(lap_row i, lap_col j) { + return data_[static_cast(i) * dim8_ + j]; + } + + [[nodiscard]] const cost& operator()(lap_row i, lap_col j) const { + return data_[static_cast(i) * dim8_ + j]; + } + + [[nodiscard]] const cost& row0(lap_col j) const { + return data_[j]; + } + + [[nodiscard]] const cost& entry0(lap_row i) const { + return data_[static_cast(i) * dim8_]; + } + + [[nodiscard]] cost* row(lap_row i) { + return &data_[static_cast(i) * dim8_]; + } + + [[nodiscard]] const cost* row(lap_row i) const { + return &data_[static_cast(i) * dim8_]; + } + + [[nodiscard]] cost* col(lap_col i) { + return &t_data_[static_cast(i) * dim8_]; + } + + [[nodiscard]] const cost* col(lap_col i) const { + return &t_data_[static_cast(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(i) * dim8_ + j] = value; + t_data_[static_cast(j) * dim8_ + i] = value; + } + + void markTransposed() noexcept { transposed_ = true; } + + // ---- Transpose ---- + + void makeTranspose() noexcept { + if (transposed_) return; + + const cost* __restrict__ data_ptr = data_.data(); + cost* __restrict__ t_data_ptr = t_data_.data(); + +#if defined(__GNUC__) || defined(__clang__) + data_ptr = + static_cast(__builtin_assume_aligned(data_ptr, 64)); + t_data_ptr = + static_cast(__builtin_assume_aligned(t_data_ptr, 64)); +#endif + for (std::size_t i = 0; i < dim_; i += BLOCK_SIZE) { + for (std::size_t j = 0; j < dim_; j += BLOCK_SIZE) { + for (std::size_t r = i; + r < std::min(i + BLOCK_SIZE, dim_); ++r) { + for (std::size_t c = j; + c < std::min(j + BLOCK_SIZE, dim_); ++c) { + t_data_ptr[c * dim8_ + r] = data_ptr[r * dim8_ + c]; + } + } + } + } + transposed_ = true; + } + + void makeUntranspose() noexcept { + const cost* __restrict__ data_ptr = t_data_.data(); + cost* __restrict__ t_data_ptr = data_.data(); + +#if defined(__GNUC__) || defined(__clang__) + data_ptr = + static_cast(__builtin_assume_aligned(data_ptr, 64)); + t_data_ptr = + static_cast(__builtin_assume_aligned(t_data_ptr, 64)); +#endif + for (std::size_t i = 0; i < dim_; i += BLOCK_SIZE) { + for (std::size_t j = 0; j < dim_; j += BLOCK_SIZE) { + for (std::size_t r = i; + r < std::min(i + BLOCK_SIZE, dim_); ++r) { + for (std::size_t c = j; + c < std::min(j + BLOCK_SIZE, dim_); ++c) { + t_data_ptr[c * dim8_ + r] = data_ptr[r * dim8_ + c]; + } + } + } + } + } + + // ---- Padding ---- + + void padTrAfterCol(lap_row start_row, cost value) { + std::size_t start_index = static_cast(start_row) * dim8_; + std::size_t end_index = dim_ * dim8_; + std::fill(t_data_.begin() + start_index, + t_data_.begin() + end_index, value); + } + + void padAfterRow(lap_row start_row, cost value) { + std::size_t start_index = static_cast(start_row) * dim8_; + std::size_t end_index = dim_ * dim8_; + std::fill(data_.begin() + start_index, + data_.begin() + end_index, value); + } + + void padTrColAfterRow(const lap_row r, const lap_col start_col, + const cost value) { + std::size_t r_offset = r * dim8_; + std::size_t actual_start_col = static_cast(start_col); + std::size_t start_index = r_offset + actual_start_col; + std::size_t end_index = start_index + dim_ - actual_start_col; + std::fill(t_data_.begin() + start_index, + t_data_.begin() + end_index, value); + } + + void padRowAfterCol(const lap_row r, const lap_col start_col, + const cost value) { + std::size_t r_offset = r * dim8_; + std::size_t actual_start_col = static_cast(start_col); + std::size_t start_index = r_offset + actual_start_col; + std::size_t end_index = start_index + dim_ - actual_start_col; + std::fill(data_.begin() + start_index, + data_.begin() + end_index, value); + } + + // ---- Search ---- + + std::pair findColMin(lap_col j, + lap_row search_dim = -1) { + makeTranspose(); + const cost* col_data = col(j); + const lap_row n = + (search_dim < 0) ? static_cast(dim_) : search_dim; + const auto min_ptr = std::min_element(col_data, col_data + n); + return {*min_ptr, + static_cast(std::distance(col_data, min_ptr))}; + } + + std::tuple findRowSubmin( + const lap_row* i, const std::vector& v + ) const { + assert(dim_ > 1); + + const cost* __restrict__ row_i = row(*i); + const lap_col dim = static_cast(dim_); + const cost* __restrict__ v_ptr = v.data(); + + const cost h0 = row_i[0] - v_ptr[0]; + const cost h1 = row_i[1] - v_ptr[1]; + + cost min_val, submin_val; + lap_col min_idx, submin_idx; + + if (h0 > h1) { + min_val = h1; submin_val = h0; + min_idx = 1; submin_idx = 0; + } else { + min_val = h0; submin_val = h1; + min_idx = 0; submin_idx = 1; + } + + const lap_col j_limit = + (dim < 4 ? 0 : static_cast(dim - 3)); + + for (lap_col j = 2; j < j_limit; j += 4) { + assert(BLOCK_SIZE >= 4); + const cost h0 = row_i[j] - v_ptr[j]; + const cost h1 = row_i[j + 1] - v_ptr[j + 1]; + const cost h2 = row_i[j + 2] - v_ptr[j + 2]; + const cost h3 = row_i[j + 3] - v_ptr[j + 3]; + if (h0 < submin_val) { + if (h0 < min_val) { + submin_val = min_val; min_val = h0; + submin_idx = min_idx; min_idx = j; + } else { submin_val = h0; submin_idx = j; } + } + if (h1 < submin_val) { + if (h1 < min_val) { + submin_val = min_val; min_val = h1; + submin_idx = min_idx; min_idx = j + 1; + } else { submin_val = h1; submin_idx = j + 1; } + } + if (h2 < submin_val) { + if (h2 < min_val) { + submin_val = min_val; min_val = h2; + submin_idx = min_idx; min_idx = j + 2; + } else { submin_val = h2; submin_idx = j + 2; } + } + if (h3 < submin_val) { + if (h3 < min_val) { + submin_val = min_val; min_val = h3; + submin_idx = min_idx; min_idx = j + 3; + } else { submin_val = h3; submin_idx = j + 3; } + } + } + for (lap_col j = j_limit; j < dim; ++j) { + const cost h = row_i[j] - v_ptr[j]; + if (h < submin_val) { + if (h < min_val) { + submin_val = min_val; min_val = h; + submin_idx = min_idx; min_idx = j; + } else { submin_val = h; submin_idx = j; } + } + } + return {min_val, submin_val, min_idx, submin_idx}; + } + +}; + +} // namespace TreeDist + +#endif // TREEDIST_COST_MATRIX_H_ diff --git a/dev/rogue_bound/ana/lap.h b/dev/rogue_bound/ana/lap.h new file mode 100644 index 00000000..673cfa22 --- /dev/null +++ b/dev/rogue_bound/ana/lap.h @@ -0,0 +1,35 @@ +#ifndef TREEDIST_LAP_H_ +#define TREEDIST_LAP_H_ + +// LAP (Linear Assignment Problem) — Jonker–Volgenant declarations. +// +// Provides the lap() function signature. The implementation lives in +// lap_impl.h, guarded by TREEDIST_LAP_IMPLEMENTATION. Downstream +// packages should include lap_impl.h in exactly one translation unit +// with the guard defined. + +#include "lap_scratch.h" + +namespace TreeDist { + + // Primary overload: caller supplies pre-allocated scratch. + extern cost lap(lap_row dim, + CostMatrix& input_cost, + std::vector& rowsol, + std::vector& colsol, + bool allow_interrupt, + LapScratch& scratch); + + // Convenience overload: creates a temporary scratch. + inline cost lap(lap_row dim, + CostMatrix& input_cost, + std::vector& rowsol, + std::vector& colsol, + bool allow_interrupt = true) { + LapScratch scratch; + return lap(dim, input_cost, rowsol, colsol, allow_interrupt, scratch); + } + +} // namespace TreeDist + +#endif // TREEDIST_LAP_H_ diff --git a/dev/rogue_bound/ana/lap_hooked.h b/dev/rogue_bound/ana/lap_hooked.h new file mode 100644 index 00000000..7a07b129 --- /dev/null +++ b/dev/rogue_bound/ana/lap_hooked.h @@ -0,0 +1,239 @@ +// Copy of TreeDist::lap() (inst/include/TreeDist/lap_impl.h) with +// - a hook called after each phase (with the current column duals v), and +// - an optional warm start from given column duals. +// Phases: 0 = column reduction, 1 = reduction transfer, 2 = ARR, +// 3 = after each augmentation (arg = augmentations done, of total). +#pragma once +#include "lap.h" + +namespace ana { +using namespace TreeDist; + +inline bool nontrivially_less_than(cost a, cost b) noexcept { + return a + ((a > ROUND_PRECISION) ? 8 : 0) < b; +} + +struct LapWork { + long long n_free_after_init = 0; // free rows entering ARR + long long n_free_after_arr = 0; // free rows entering augmentation + long long arr_steps = 0; // findRowSubmin calls + long long aug_inner = 0; // Dijkstra inner-loop column visits + long long aug_scan = 0; // Dijkstra min-scan column visits + unsigned long long t_init = 0, t_arr = 0, t_aug = 0; +}; + +template +cost lap_hooked(const lap_row dim, CostMatrix &input_cost, + std::vector &rowsol, std::vector &colsol, + LapScratch &scratch, const cost *v_init, Hook &&hook, + LapWork &work, const int arr_loops = 2) { + unsigned long long tick = __rdtsc(); + lap_row num_free = 0; + scratch.ensure(dim); + auto &v = scratch.v; + auto &matches = scratch.matches; + std::fill(matches.begin(), matches.begin() + dim, 0); + const cost *__restrict__ v_ptr = v.data(); + auto &freeunassigned = scratch.freeunassigned; + + if (!v_init) { + // COLUMN REDUCTION + for (lap_col j = dim; j--;) { + const auto [min, imin] = input_cost.findColMin(j); + v[j] = min; + ++matches[imin]; + if (matches[imin] == 1) { + rowsol[imin] = j; + colsol[j] = imin; + } else if (v_ptr[j] < v_ptr[rowsol[imin]]) { + const lap_col j1 = rowsol[imin]; + rowsol[imin] = j; + colsol[j] = imin; + colsol[j1] = -1; + } else { + colsol[j] = -1; + } + } + hook(0, v, 0, 0); + + // REDUCTION TRANSFER + for (lap_row i = 0; i < dim; ++i) { + if (matches[i] == 0) { + freeunassigned[num_free++] = i; + } else if (matches[i] == 1) { + const lap_col j1 = rowsol[i]; + const cost *row_i = input_cost.row(i); + cost min_cost; + if (j1 == 0) { + min_cost = row_i[1] - v_ptr[1]; + for (lap_col j = 2; j < dim; ++j) { + const cost reduced_cost = row_i[j] - v_ptr[j]; + if (reduced_cost < min_cost) min_cost = reduced_cost; + } + } else { + min_cost = row_i[0] - v_ptr[0]; + for (lap_col j = 1; j < dim; ++j) { + if (j == j1) continue; + const cost reduced_cost = row_i[j] - v_ptr[j]; + if (reduced_cost < min_cost) min_cost = reduced_cost; + } + } + v[j1] -= min_cost; + } + } + hook(1, v, 0, 0); + } else { + // WARM START: v given; assign each row to its cheapest reduced column + // if that column is still free, so assigned edges are tight. + for (lap_col j = 0; j < dim; ++j) { + v[j] = v_init[j]; + colsol[j] = -1; + } + for (lap_row i = 0; i < dim; ++i) { + const cost *row_i = input_cost.row(i); + cost best = row_i[0] - v_ptr[0]; + lap_col best_j = 0; + for (lap_col j = 1; j < dim; ++j) { + const cost h = row_i[j] - v_ptr[j]; + if (h < best) { best = h; best_j = j; } + } + if (colsol[best_j] < 0) { + colsol[best_j] = i; + rowsol[i] = best_j; + } else { + freeunassigned[num_free++] = i; + } + } + hook(1, v, 0, 0); + } + work.n_free_after_init += num_free; + { unsigned long long t = __rdtsc(); work.t_init += t - tick; tick = t; } + + // AUGMENTING ROW REDUCTION + int loopcnt = 0; + while (loopcnt < arr_loops) { + ++loopcnt; + lap_row previous_num_free = num_free; + num_free = 0; + lap_row k = 0; + while (k < previous_num_free) { + const lap_row i = freeunassigned[k++]; + ++work.arr_steps; + const auto [umin, usubmin, min_idx, j2] = input_cost.findRowSubmin(&i, v); + lap_col j1 = min_idx; + lap_row i0 = colsol[j1]; + const bool strictly_less = nontrivially_less_than(umin, usubmin); + if (strictly_less) { + v[j1] -= (usubmin - umin); + } else if (i0 > -1) { + j1 = j2; + i0 = colsol[j2]; + } + rowsol[i] = j1; + colsol[j1] = i; + if (i0 > -1) { + if (strictly_less) { + freeunassigned[--k] = i0; + } else { + freeunassigned[num_free++] = i0; + } + } + } + } + hook(2, v, 0, num_free); + work.n_free_after_arr += num_free; + { unsigned long long t = __rdtsc(); work.t_arr += t - tick; tick = t; } + + cost *__restrict__ d_ptr = scratch.d.data(); + lap_row *__restrict__ pred_ptr = scratch.predecessor.data(); + lap_col *__restrict__ cl_ptr = scratch.col_list.data(); + + for (lap_row f = 0; f < num_free; ++f) { + bool unassignedfound = false; + lap_row free_row = freeunassigned[f]; + const cost *__restrict__ free_row_cost = input_cost.row(free_row); + lap_col endofpath = 0; + lap_col last = 0; + lap_row i; + lap_col j1; + for (lap_col j = 0; j < dim; ++j) { + d_ptr[j] = free_row_cost[j] - v_ptr[j]; + pred_ptr[j] = free_row; + cl_ptr[j] = j; + } + cost min = 0; + lap_col low = 0; + lap_col up = 0; + do { + if (up == low) { + last = low - 1; + min = d_ptr[cl_ptr[up++]]; + for (lap_dim k = up; k < dim; ++k) { + ++work.aug_scan; + const lap_col j = cl_ptr[k]; + const cost h = d_ptr[j]; + if (h <= min) { + if (h < min) { + up = low; + min = h; + } + cl_ptr[k] = cl_ptr[up]; + cl_ptr[up++] = j; + } + } + for (lap_dim k = low; k < up; ++k) { + if (colsol[cl_ptr[k]] < 0) { + endofpath = cl_ptr[k]; + unassignedfound = true; + break; + } + } + } + if (!unassignedfound) { + j1 = cl_ptr[low++]; + i = colsol[j1]; + const cost *__restrict__ row_i = input_cost.row(i); + const cost h = row_i[j1] - v_ptr[j1] - min; + for (lap_dim k = up; k < dim; ++k) { + ++work.aug_inner; + const lap_col j = cl_ptr[k]; + cost v2 = row_i[j] - v_ptr[j] - h; + if (v2 < d_ptr[j]) { + pred_ptr[j] = i; + if (v2 == min) { + if (colsol[j] < 0) { + endofpath = j; + unassignedfound = true; + break; + } else { + cl_ptr[k] = cl_ptr[up]; + cl_ptr[up++] = j; + } + } + d_ptr[j] = v2; + } + } + } + } while (!unassignedfound); + for (lap_dim k = 0; k <= last; ++k) { + j1 = cl_ptr[k]; + v[j1] += d_ptr[j1] - min; + } + do { + i = pred_ptr[endofpath]; + colsol[endofpath] = i; + j1 = endofpath; + endofpath = rowsol[i]; + rowsol[i] = j1; + } while (i != free_row); + hook(3, v, f + 1, num_free); + } + + { unsigned long long t = __rdtsc(); work.t_aug += t - tick; tick = t; } + cost lapcost = 0; + for (lap_dim i = 0; i < dim; ++i) { + lapcost += input_cost(i, rowsol[i]); + } + return lapcost; +} +} // namespace ana diff --git a/dev/rogue_bound/ana/lap_scratch.h b/dev/rogue_bound/ana/lap_scratch.h new file mode 100644 index 00000000..4a9ca43b --- /dev/null +++ b/dev/rogue_bound/ana/lap_scratch.h @@ -0,0 +1,44 @@ +#ifndef TREEDIST_LAP_SCRATCH_H_ +#define TREEDIST_LAP_SCRATCH_H_ + +// Reusable heap storage for one thread's LAPJV workspace. +// +// Pass to the scratch-taking overload of lap() to eliminate per-call +// std::vector allocations. In a serial context construct once before the +// loop; in an OpenMP context allocate one per thread and index by +// omp_get_thread_num(). All vectors are grown lazily; never shrunk. + +#include "cost_matrix.h" + +namespace TreeDist { + +struct LapScratch { + std::vector v; + std::vector matches; + std::vector freeunassigned; + std::vector col_list; + std::vector d; + std::vector predecessor; + std::vector rowsol; + std::vector colsol; + CostMatrix score_pool; + CostMatrix small_pool; + + void ensure(int dim) noexcept { + const int padded = static_cast( + ((static_cast(dim) + BLOCK_SIZE - 1) + / BLOCK_SIZE) * BLOCK_SIZE); + if (static_cast(v.size()) < padded) v.resize(padded); + if (static_cast(matches.size()) < dim) matches.resize(dim); + if (static_cast(freeunassigned.size()) < dim) freeunassigned.resize(dim); + if (static_cast(col_list.size()) < dim) col_list.resize(dim); + if (static_cast(d.size()) < dim) d.resize(dim); + if (static_cast(predecessor.size()) < dim) predecessor.resize(dim); + if (static_cast(rowsol.size()) < dim) rowsol.resize(dim); + if (static_cast(colsol.size()) < dim) colsol.resize(dim); + } +}; + +} // namespace TreeDist + +#endif // TREEDIST_LAP_SCRATCH_H_ diff --git a/dev/rogue_bound/ana/stub/TreeTools/SplitList.h b/dev/rogue_bound/ana/stub/TreeTools/SplitList.h new file mode 100644 index 00000000..947e88ab --- /dev/null +++ b/dev/rogue_bound/ana/stub/TreeTools/SplitList.h @@ -0,0 +1,9 @@ +// Stub for analysis builds: the constants TreeDist/types.h needs, no Rcpp. +#pragma once +#include +using splitbit = uint_fast64_t; +#define R_BIN_SIZE int16_t(8) +#define SL_BIN_SIZE int16_t(64) +#define SL_MAX_BINS int16_t(511) +#define SL_MAX_TIPS (SL_BIN_SIZE * SL_MAX_BINS) +#define SL_MAX_SPLITS (SL_MAX_TIPS - 3) diff --git a/dev/rogue_bound/ana/types.h b/dev/rogue_bound/ana/types.h new file mode 100644 index 00000000..9023ecf5 --- /dev/null +++ b/dev/rogue_bound/ana/types.h @@ -0,0 +1,36 @@ +#ifndef TREEDIST_TYPES_H_ +#define TREEDIST_TYPES_H_ + +// TreeDist public type definitions. +// +// Rcpp-free: uses only standard headers plus TreeTools (for SL_MAX_SPLITS). +// Downstream packages include via: #include + +#include +#include +#include "stub/TreeTools/SplitList.h" + +namespace TreeDist { + + using int16 = int_fast16_t; + using int32 = int_fast32_t; + using cost = int_fast64_t; + + // Primary type for split/tip/bin counters. int32 benchmarks faster than int16 + // (avoids sign-extension when mixing with SplitList int32 members). + using split_int = int32; + + using lap_dim = int; + using lap_row = lap_dim; + using lap_col = lap_dim; + + constexpr cost BIG = + (std::numeric_limits::max)() / SL_MAX_SPLITS; + + constexpr cost ROUND_PRECISION = 2048 * 2048; + + constexpr std::size_t BLOCK_SIZE = 16; + +} // namespace TreeDist + +#endif // TREEDIST_TYPES_H_ diff --git a/dev/rogue_bound/bench.R b/dev/rogue_bound/bench.R new file mode 100644 index 00000000..0a6d3e6a --- /dev/null +++ b/dev/rogue_bound/bench.R @@ -0,0 +1,16 @@ +# Usage: Rscript bench.R [depth] +# Prints min and median elapsed time over `reps` runs, and a checksum. +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +suppressMessages(library(TreeTools)) +wl <- readRDS("/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/workloads.rds") +pairs <- wl[[args[2]]] +reps <- as.integer(args[3]) +depth <- if (length(args) > 3) as.integer(args[4]) else 1L +options(sprDepth = depth) +Run <- function() vapply(pairs, function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +r <- Run() +el <- replicate(reps, system.time(Run())[["elapsed"]]) +cat(sprintf("%-10s %-8s d=%d min %.3f med %.3f sum=%g\n", basename(args[1]), + args[2], depth, min(el), median(el), sum(r))) diff --git a/dev/rogue_bound/bench_results.txt b/dev/rogue_bound/bench_results.txt new file mode 100644 index 00000000..a8805a20 --- /dev/null +++ b/dev/rogue_bound/bench_results.txt @@ -0,0 +1,49 @@ +lib_base rand30 d=1 min 0.464 med 0.605 sum=811 +lib_arr rand30 d=1 min 0.433 med 0.459 sum=811 +lib_new2 rand30 d=1 min 0.343 med 0.457 sum=811 +lib_base near80 d=1 min 0.049 med 0.052 sum=110 +lib_arr near80 d=1 min 0.049 med 0.053 sum=110 +lib_new2 near80 d=1 min 0.047 med 0.052 sum=110 +lib_base near150 d=1 min 0.053 med 0.058 sum=80 +lib_arr near150 d=1 min 0.049 med 0.052 sum=80 +lib_new2 near150 d=1 min 0.055 med 0.056 sum=80 +lib_base mid80 d=1 min 0.242 med 0.264 sum=166 +lib_arr mid80 d=1 min 0.222 med 0.233 sum=166 +lib_new2 mid80 d=1 min 0.219 med 0.245 sum=166 +lib_base mid150 d=1 min 1.356 med 1.442 sum=161 +lib_arr mid150 d=1 min 1.300 med 1.408 sum=161 +lib_new2 mid150 d=1 min 1.001 med 1.101 sum=161 +lib_base rand80 d=1 min 4.089 med 4.280 sum=777 +lib_arr rand80 d=1 min 3.565 med 3.669 sum=777 +lib_new2 rand80 d=1 min 2.162 med 2.308 sum=777 +lib_base rand150 d=1 min 18.244 med 18.905 sum=520 +lib_arr rand150 d=1 min 16.858 med 17.005 sum=520 +lib_new2 rand150 d=1 min 6.396 med 6.539 sum=520 +lib_base rand30 d=2 min 3.698 med 3.998 sum=800 +lib_arr rand30 d=2 min 3.350 med 3.407 sum=800 +lib_new2 rand30 d=2 min 3.131 med 3.161 sum=800 +lib_base rand30 d=1 min 0.423 med 0.489 sum=811 +lib_arr rand30 d=1 min 0.380 med 0.437 sum=811 +lib_new2 rand30 d=1 min 0.343 med 0.378 sum=811 +lib_base near80 d=1 min 0.048 med 0.053 sum=110 +lib_arr near80 d=1 min 0.045 med 0.050 sum=110 +lib_new2 near80 d=1 min 0.047 med 0.053 sum=110 +lib_base near150 d=1 min 0.054 med 0.055 sum=80 +lib_arr near150 d=1 min 0.054 med 0.056 sum=80 +lib_new2 near150 d=1 min 0.056 med 0.058 sum=80 +lib_base mid80 d=1 min 0.243 med 0.249 sum=166 +lib_arr mid80 d=1 min 0.236 med 0.241 sum=166 +lib_new2 mid80 d=1 min 0.226 med 0.275 sum=166 +lib_base mid150 d=1 min 1.307 med 1.616 sum=161 +lib_arr mid150 d=1 min 1.310 med 1.402 sum=161 +lib_new2 mid150 d=1 min 1.062 med 1.126 sum=161 +lib_base rand80 d=1 min 3.997 med 4.170 sum=777 +lib_arr rand80 d=1 min 3.442 med 3.566 sum=777 +lib_new2 rand80 d=1 min 1.944 med 2.214 sum=777 +lib_base rand150 d=1 min 18.209 med 18.458 sum=520 +lib_arr rand150 d=1 min 15.795 med 16.962 sum=520 +lib_new2 rand150 d=1 min 6.682 med 7.497 sum=520 +lib_base rand30 d=2 min 3.659 med 3.667 sum=800 +lib_arr rand30 d=2 min 3.293 med 3.548 sum=800 +lib_new2 rand30 d=2 min 3.218 med 3.233 sum=800 +DONE diff --git a/dev/rogue_bound/bench_summary.R b/dev/rogue_bound/bench_summary.R new file mode 100644 index 00000000..98432174 --- /dev/null +++ b/dev/rogue_bound/bench_summary.R @@ -0,0 +1,16 @@ +# Min time per (lib, workload, depth) across interleaved rounds; % vs base. +x <- read.table(text = grep("^lib_", readLines( + "/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/bench_results.txt"), + value = TRUE), stringsAsFactors = FALSE) +names(x) <- c("lib", "wl", "depth", "minLab", "min", "medLab", "med", "sum") +x$key <- paste(x$wl, x$depth) +mins <- tapply(x$min, list(x$key, x$lib), min) +runs <- tapply(x$min, list(x$key, x$lib), length) +out <- data.frame(base = mins[, "lib_base"], arr = mins[, "lib_arr"], + new2 = mins[, "lib_new2"], + arrPct = round(100 * (mins[, "lib_arr"] / mins[, "lib_base"] - 1)), + new2Pct = round(100 * (mins[, "lib_new2"] / mins[, "lib_base"] - 1)), + rounds = runs[, "lib_base"]) +print(out) +sums <- tapply(x$sum, x$key, function(s) length(unique(s))) +cat("checksums agree across libs:", all(sums == 1), "\n") diff --git a/dev/rogue_bound/compare.R b/dev/rogue_bound/compare.R new file mode 100644 index 00000000..a4184c62 --- /dev/null +++ b/dev/rogue_bound/compare.R @@ -0,0 +1,52 @@ +# Usage: Rscript compare.R +# Checks, for each spr_rogue_scores() call, that both versions agree on the +# minimum score and on the set of candidates attaining it, and that every +# candidate scored exactly by either version (score within tolerance of the +# minimum) has the same score in both. +args <- commandArgs(TRUE) +a <- readRDS(args[1]) +b <- readRDS(args[2]) +stopifnot(length(a) == length(b)) +Check <- function(x, y, label) { + if (length(x) != length(y)) return(sprintf("%s: length", label)) + if (!length(x)) return(NULL) + fin <- is.finite(x) & is.finite(y) + if (!identical(is.finite(x), is.finite(y))) return(sprintf("%s: finiteness", label)) + if (!any(fin)) return(NULL) + mx <- min(x[fin]); my <- min(y[fin]) + tol <- 1e-9 * max(1, abs(mx)) + if (abs(mx - my) > tol) return(sprintf("%s: min %.15g vs %.15g", label, mx, my)) + # Candidates tied with the minimum (exact equality, as SPRDist uses) + if (!identical(which(x == mx), which(y == my))) { + near <- abs(x - mx) <= tol | abs(y - my) <= tol + return(sprintf("%s: argmin set differs (%s vs %s); near-ties: %d", + label, paste(which(x == mx), collapse = " "), + paste(which(y == my), collapse = " "), sum(near))) + } + # Every candidate close to the minimum must be exact in both + close <- fin & (x <= mx + 0.5 | y <= my + 0.5) + both <- close & x <= mx + 1e-6 & y <= my + 1e-6 + if (any(abs(x[both] - y[both]) > tol)) return(sprintf("%s: near-min scores differ", label)) + # A candidate rejected by one version must be rejected (> min) by the other + if (any(x[fin] <= mx & y[fin] > my) || any(y[fin] <= my & x[fin] > mx)) { + return(sprintf("%s: rejection mismatch", label)) + } + NULL +} +problems <- character(0) +nCand <- 0 +for (i in seq_along(a)) { + p1 <- Check(a[[i]]$scores, b[[i]]$scores, sprintf("call %d singles", i)) + p2 <- Check(a[[i]]$pair_scores, b[[i]]$pair_scores, sprintf("call %d pairs", i)) + nCand <- nCand + length(a[[i]]$scores) + length(a[[i]]$pair_scores) + problems <- c(problems, p1, p2) +} +cat("calls:", length(a), " candidates:", nCand, " problems:", length(problems), "\n") +if (length(problems)) print(head(problems, 20)) +# How often exact scores differ at all (machine-precision level) +diffs <- unlist(lapply(seq_along(a), function(i) { + x <- a[[i]]$scores; y <- b[[i]]$scores + m <- min(x); keep <- is.finite(x) & x == m & y == min(y) + abs(x[keep] - y[keep]) +})) +cat("max |diff| among minimum scores:", max(c(0, diffs)), "\n") diff --git a/dev/rogue_bound/instrumentation_dump.patch b/dev/rogue_bound/instrumentation_dump.patch new file mode 100644 index 00000000..bbb3fd87 --- /dev/null +++ b/dev/rogue_bound/instrumentation_dump.patch @@ -0,0 +1,339 @@ +diff --git a/src/instr.h b/src/instr.h +new file mode 100644 +index 0000000..558e2ab +--- /dev/null ++++ b/src/instr.h +@@ -0,0 +1,32 @@ ++// Instrumentation for bound research (scratch only). ++#ifndef TREEDIST_INSTR_H_ ++#define TREEDIST_INSTR_H_ ++#include ++#include ++#include ++ ++namespace instr { ++enum Timer { ++ T_TOTAL, T_KAR, T_SPLITS, T_FMS, T_BCID, T_MCI_SETUP, T_MCI_FILL, ++ T_MCI_COPY, T_MCI_LAP, T_MCI_TAIL, T_HALF_KAR, T_PARENT, N_TIMER ++}; ++enum Count { ++ C_CALLS, C_ENTROPY, C_ROW, C_EXACT, C_LAP, C_LAP_BELOW, C_WHOLE, C_HALF1, ++ C_HALF2, C_CANDIDATES, C_EXTRA_EXIT, N_COUNT ++}; ++inline uint64_t ticks[N_TIMER] = {}; ++inline uint64_t counts[N_COUNT] = {}; ++inline FILE *dump = nullptr; ++inline int cand_id = -1; ++inline int call_kind = -1; ++inline bool lap_below_last = false; ++ ++struct Scope { ++ int t; ++ uint64_t s; ++ explicit Scope(int timer) : t(timer), s(__rdtsc()) {} ++ ~Scope() { ticks[t] += __rdtsc() - s; } ++}; ++} ++ ++#endif +diff --git a/src/spr.cpp b/src/spr.cpp +index 7651639..ac25132 100644 +--- a/src/spr.cpp ++++ b/src/spr.cpp +@@ -5,6 +5,8 @@ + #include // for root_on_node + #include // for SplitList, count_bits + #include "reduce_tree.h" ++#include "instr.h" ++#include + #include // for log2 + #include // for sort, lower_bound, equal + #include +@@ -251,6 +253,20 @@ extern double mutual_clustering_bounded_score(const RawMatrix &x, + const double min_score); + + namespace { ++std::unordered_map g_parent_index; ++const CharacterVector *g_cur_labels = nullptr; ++ ++void DumpSplits(const RawMatrix &m) { ++ const int nr = m.nrow(), nc = m.ncol(); ++ int hdr[2] = {nr, nc}; ++ fwrite(hdr, sizeof(int), 2, instr::dump); ++ for (int r = 0; r != nr; ++r) { ++ for (int c = 0; c != nc; ++c) { ++ const uint8_t b = m(r, c); ++ fwrite(&b, 1, 1, instr::dump); ++ } ++ } ++} + + // .Machine$double.eps ^ 0.5 + constexpr double SQRT_EPS = 1.4901161193847656e-08; +@@ -330,11 +346,31 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { + // numerical noise, may instead return a lower bound that exceeds `cutoff`. + double BoundedCID(const TreeSplits &splits1, const TreeSplits &splits2, + const int n_tip, const double cutoff) { ++ instr::Scope scope(instr::T_BCID); + const double info = splits1.entropy + splits2.entropy; + const double tol = std::max(1.0, info) * SQRT_EPS; + const double mci = n_tip < 4 ? 0 : + mutual_clustering_bounded_score(splits1.bits, splits2.bits, n_tip, + (info - cutoff) / 2 - tol); ++ if (instr::dump && n_tip >= 4) { ++ const double min_score = (info - cutoff) / 2 - tol; ++ const double exact = mutual_clustering_bounded_score( ++ splits1.bits, splits2.bits, n_tip, R_NegInf); ++ int hdr[4] = {2, instr::cand_id, instr::call_kind, n_tip}; ++ fwrite(hdr, sizeof(int), 4, instr::dump); ++ std::vector map(n_tip, -1); ++ if (g_cur_labels) { ++ for (int t = 0; t != n_tip; ++t) { ++ auto it = g_parent_index.find(STRING_ELT(*g_cur_labels, t)); ++ if (it != g_parent_index.end()) map[t] = it->second; ++ } ++ } ++ fwrite(map.data(), sizeof(int), n_tip, instr::dump); ++ double vals[4] = {min_score, exact, mci, info}; ++ fwrite(vals, sizeof(double), 4, instr::dump); ++ DumpSplits(splits1.bits); ++ DumpSplits(splits2.bits); ++ } + const double ret = info - mci - mci; + return ret < tol ? 0 : ret; + } +@@ -342,8 +378,16 @@ double BoundedCID(const TreeSplits &splits1, const TreeSplits &splits2, + double BoundedCID(const List &tree1, const List &tree2, const double cutoff) { + const CharacterVector labels = tree1["tip.label"]; + const int n_tip = labels.size(); +- return BoundedCID(GetSplits(tree1["edge"], n_tip), +- GetSplits(tree2["edge"], n_tip), n_tip, cutoff); ++ g_cur_labels = &labels; ++ TreeSplits s1, s2; ++ { ++ instr::Scope scope(instr::T_SPLITS); ++ s1 = GetSplits(tree1["edge"], n_tip); ++ s2 = GetSplits(tree2["edge"], n_tip); ++ } ++ const double ret = BoundedCID(s1, s2, n_tip, cutoff); ++ g_cur_labels = nullptr; ++ return ret; + } + + // Index of the first split in `x` that also occurs in `table`, or -1; +@@ -426,7 +470,10 @@ List KeepAndReduce(const List &tree1, const List &tree2, + double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + const CharacterVector &labels, + const LogicalVector &keep, const double cutoff) { ++ ++instr::counts[instr::C_CANDIDATES]; ++ uint64_t t0 = __rdtsc(); + const List outcome = KeepAndReduce(pre1, pre2, labels, keep); ++ { const uint64_t t = __rdtsc(); instr::ticks[instr::T_KAR] += t - t0; t0 = t; } + if (Rf_isNull(outcome[0])) { + return R_NegInf; + } +@@ -436,10 +483,17 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + const int n_tip = out_labels.size(); + const TreeSplits splits1 = GetSplits(out1["edge"], n_tip); + const TreeSplits splits2 = GetSplits(out2["edge"], n_tip); ++ { const uint64_t t = __rdtsc(); instr::ticks[instr::T_SPLITS] += t - t0; t0 = t; } + const int first_match = FirstMatchingSplit(splits1.bits, splits2.bits, + n_tip); ++ { const uint64_t t = __rdtsc(); instr::ticks[instr::T_FMS] += t - t0; t0 = t; } + if (first_match < 0) { +- return BoundedCID(splits1, splits2, n_tip, cutoff); ++ ++instr::counts[instr::C_WHOLE]; ++ instr::call_kind = 0; ++ g_cur_labels = &out_labels; ++ const double ret = BoundedCID(splits1, splits2, n_tip, cutoff); ++ g_cur_labels = nullptr; ++ return ret; + } + + // Divide at the shared split, anchoring each half to the other +@@ -456,18 +510,26 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + } + } + } ++ t0 = __rdtsc(); + const List sub1 = KeepAndReduce(out1, out2, half1); ++ { const uint64_t t = __rdtsc(); instr::ticks[instr::T_HALF_KAR] += t - t0; t0 = t; } + if (Rf_isNull(sub1[0])) { + Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov + } ++ ++instr::counts[instr::C_HALF1]; ++ instr::call_kind = 1; + const double score1 = BoundedCID(sub1[0], sub1[1], cutoff); + if (score1 > cutoff) { + return score1; + } ++ t0 = __rdtsc(); + const List sub2 = KeepAndReduce(out1, out2, half2); ++ { const uint64_t t = __rdtsc(); instr::ticks[instr::T_HALF_KAR] += t - t0; t0 = t; } + if (Rf_isNull(sub2[0])) { + Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov + } ++ ++instr::counts[instr::C_HALF2]; ++ instr::call_kind = 2; + return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1); + } + +@@ -480,8 +542,21 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + // rather than their exact score. + // [[Rcpp::export]] + List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { ++ instr::Scope total_scope(instr::T_TOTAL); + const CharacterVector labels = tree1["tip.label"]; + const int n_tip = labels.size(); ++ if (instr::dump) { ++ g_parent_index.clear(); ++ for (int t = 0; t != n_tip; ++t) g_parent_index[STRING_ELT(labels, t)] = t; ++ const TreeSplits p1 = GetSplits(tree1["edge"], n_tip); ++ const TreeSplits p2 = GetSplits(tree2["edge"], n_tip); ++ int hdr[2] = {1, n_tip}; ++ fwrite(hdr, sizeof(int), 2, instr::dump); ++ double info = p1.entropy + p2.entropy; ++ fwrite(&info, sizeof(double), 1, instr::dump); ++ DumpSplits(p1.bits); ++ DumpSplits(p2.bits); ++ } + const IntegerMatrix pre1 = reverse(tree1["edge"]); + const IntegerMatrix pre2 = reverse(tree2["edge"]); + NumericVector scores(n_tip); +@@ -491,6 +566,7 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { + for (int i = 0; i != n_tip; ++i) { + Rcpp::checkUserInterrupt(); + keep[i] = false; ++ instr::cand_id = i; + scores[i] = ScoreWithout(pre1, pre2, labels, keep, best); + keep[i] = true; + if (!R_FINITE(scores[i])) { +@@ -512,6 +588,7 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { + keep[i] = false; + for (int j = i + 1; j != n_tip; ++j) { + keep[j] = false; ++ instr::cand_id = n_tip + pair; + pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best); + keep[j] = true; + if (!R_FINITE(pair_scores[pair])) { +@@ -562,3 +639,21 @@ List spr_keep_and_reduce(const List &tree1, const List &tree2, + const LogicalVector &keep) { + return KeepAndReduce(tree1, tree2, keep); + } ++ ++// [[Rcpp::export]] ++NumericVector instr_get(const bool reset) { ++ NumericVector ret(instr::N_TIMER + instr::N_COUNT); ++ for (int i = 0; i != instr::N_TIMER; ++i) ret[i] = double(instr::ticks[i]); ++ for (int i = 0; i != instr::N_COUNT; ++i) ret[instr::N_TIMER + i] = double(instr::counts[i]); ++ if (reset) { ++ for (auto &t : instr::ticks) t = 0; ++ for (auto &c : instr::counts) c = 0; ++ } ++ return ret; ++} ++ ++// [[Rcpp::export]] ++void instr_dump(const std::string &path) { ++ if (instr::dump) { fclose(instr::dump); instr::dump = nullptr; } ++ if (!path.empty()) instr::dump = fopen(path.c_str(), "wb"); ++} +diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp +index 0dc6eae..624e8a5 100644 +--- a/src/tree_distances.cpp ++++ b/src/tree_distances.cpp +@@ -5,6 +5,7 @@ + #include /* for unique_ptr, make_unique */ + #include + #include "tree_distances.h" ++#include "instr.h" + + using Rcpp::_; + using Rcpp::IntegerVector; +@@ -417,6 +418,8 @@ static cost_matrix &BoundedPool(const int i, const std::size_t dim) { + template + List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + const int32 n_tips, const double min_score = 0) { ++ [[maybe_unused]] uint64_t tick0 = __rdtsc(); ++ if constexpr (Bounded) ++instr::counts[instr::C_CALLS]; + const SplitList a(x); + const SplitList b(y); + const bool a_has_more_splits = (a.n_splits > b.n_splits); +@@ -451,6 +454,7 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + } + upper_bound = std::min(entropy_a, entropy_b) * n_tips_reciprocal; + if (upper_bound < min_score) { ++ ++instr::counts[instr::C_ENTROPY]; + return List::create(Named("score") = upper_bound, + _["matching"] = IntegerVector(0)); + } +@@ -461,6 +465,9 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + // buffers rather than reallocating them on each call. + cost_matrix local_score(Bounded ? 0 : most_splits); + cost_matrix &score = Bounded ? BoundedPool(0, most_splits) : local_score; ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_SETUP] += t - tick0; tick0 = t; ++ } + + double exact_match_score = 0; + split_int exact_matches = 0; +@@ -530,6 +537,8 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + over_max_score - TreeDist::ic_matching(na, nA, n_tips) * + n_tips_reciprocal; + if (upper_bound < min_score) { ++ ++instr::counts[instr::C_ROW]; ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_FILL] += t - tick0; tick0 = t; + return List::create(Named("score") = upper_bound, + _["matching"] = IntegerVector(0)); + } +@@ -541,7 +550,11 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + } + } + ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_FILL] += t - tick0; tick0 = t; ++ } + if (exact_matches == b.n_splits || exact_matches == a.n_splits) { ++ if constexpr (Bounded) ++instr::counts[instr::C_EXACT]; + return List::create( + Named("score") = exact_match_score * n_tips_reciprocal, + _["matching"] = a_match); +@@ -570,10 +583,19 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + a_pos++; + } + small_score.padAfterRow(lap_dim - b_extra_splits, max_score); ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_COPY] += t - tick0; tick0 = t; ++ } + + const double lap_score = static_cast((max_score * lap_dim) - + lap(lap_dim, small_score, rowsol, colsol, !Bounded)) * over_max_score; + const double final_score = lap_score + (exact_match_score / n_tips); ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_LAP] += t - tick0; tick0 = t; ++ ++instr::counts[instr::C_LAP]; ++ instr::lap_below_last = final_score < min_score; ++ if (final_score < min_score) ++instr::counts[instr::C_LAP_BELOW]; ++ } + + std::unique_ptr lap_decode = std::make_unique(lap_dim); + split_int fuzzy_match = 0; +@@ -609,9 +631,18 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + } + } + ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_COPY] += t - tick0; tick0 = t; ++ } + const double final_score = static_cast( + (max_score * lap_dim) - lap(lap_dim, score, rowsol, colsol, !Bounded) + ) / max_score; ++ if constexpr (Bounded) { ++ const uint64_t t = __rdtsc(); instr::ticks[instr::T_MCI_LAP] += t - tick0; tick0 = t; ++ ++instr::counts[instr::C_LAP]; ++ instr::lap_below_last = final_score < min_score; ++ if (final_score < min_score) ++instr::counts[instr::C_LAP_BELOW]; ++ } + + std::vector final_matching; + final_matching.reserve(a.n_splits); diff --git a/dev/rogue_bound/make_workload.R b/dev/rogue_bound/make_workload.R new file mode 100644 index 00000000..193dea4e --- /dev/null +++ b/dev/rogue_bound/make_workload.R @@ -0,0 +1,49 @@ +# Generate fixed workloads of tree pairs for SPRDist(method = "Rogue"). +# Random pairs: two independent RandomTree(n, root = TRUE). +# Near pairs: tree2 = tree1 after `nMoves` leaf prune-and-regraft moves. +suppressMessages(library(TreeTools, lib.loc = "/root/Rlib")) + +MovedTree <- function(tree, nMoves) { + for (i in seq_len(nMoves)) { + leaf <- sample(tree[["tip.label"]], 1) + pruned <- DropTip(tree, leaf) + nEdge <- dim(pruned[["edge"]])[1] + where <- pruned[["edge"]][sample.int(nEdge, 1), 2] + tree <- AddTip(pruned, where = where, label = leaf) + } + RenumberTips(Preorder(RootTree(tree, 1)), sort(tree[["tip.label"]])) +} + +RandomPair <- function(nTip) { + t1 <- RandomTree(nTip, root = TRUE) + t2 <- RandomTree(nTip, root = TRUE) + list(t1, RenumberTips(t2, t1)) +} + +NearPair <- function(nTip, nMoves) { + t1 <- RandomTree(nTip, root = TRUE) + t1 <- RenumberTips(t1, sort(t1[["tip.label"]])) + list(t1, RenumberTips(MovedTree(t1, nMoves), t1)) +} + +set.seed(20261014) +workloads <- list( + rand30 = lapply(1:40, function(i) RandomPair(30)), + rand80 = lapply(1:12, function(i) RandomPair(80)), + rand150 = lapply(1:4, function(i) RandomPair(150)), + near80 = lapply(1:30, function(i) NearPair(80, sample(2:6, 1))), + near150 = lapply(1:16, function(i) NearPair(150, sample(3:8, 1))) +) +saveRDS(workloads, "workloads.rds") + +# Larger set for output verification (>= 200 pairs) +set.seed(4242) +verify <- c( + lapply(1:60, function(i) RandomPair(sample(8:40, 1))), + lapply(1:40, function(i) RandomPair(sample(41:100, 1))), + lapply(1:10, function(i) RandomPair(sample(101:160, 1))), + lapply(1:60, function(i) NearPair(sample(10:60, 1), sample(1:5, 1))), + lapply(1:40, function(i) NearPair(sample(61:150, 1), sample(2:8, 1))) +) +saveRDS(verify, "verify.rds") +cat(vapply(workloads, length, 1), length(verify), "\n") diff --git a/dev/rogue_bound/parent_prices_warm_start.patch b/dev/rogue_bound/parent_prices_warm_start.patch new file mode 100644 index 00000000..01d8f4ff --- /dev/null +++ b/dev/rogue_bound/parent_prices_warm_start.patch @@ -0,0 +1,764 @@ +diff --git a/src/col_prices.h b/src/col_prices.h +new file mode 100644 +index 0000000..65b4eed +--- /dev/null ++++ b/src/col_prices.h +@@ -0,0 +1,15 @@ ++#ifndef TREEDIST_COL_PRICES_H_ ++#define TREEDIST_COL_PRICES_H_ ++ ++#include ++ ++// Supplies, on demand, prices for the splits of the second tree in a bounded ++// mutual clustering calculation: one per split, in LAP cost units, each in ++// [0, BIG]. `get` may return null if no prices are available. Prices ++// affect only speed, never results. ++struct ColPrices { ++ const int_fast64_t *(*get)(void *context); ++ void *context; ++}; ++ ++#endif +diff --git a/src/lap_warm.cpp b/src/lap_warm.cpp +new file mode 100644 +index 0000000..2b61b0b +--- /dev/null ++++ b/src/lap_warm.cpp +@@ -0,0 +1,182 @@ ++// Jonker-Volgenant LAP warm-started from given column duals. ++// ++// As TreeDist::lap() (inst/include/TreeDist/lap_impl.h), except that column ++// reduction and reduction transfer are replaced by a start from the column ++// duals `v_init`: each row takes its cheapest reduced-cost column if that ++// column is still free, so every assigned edge is tight, as the augmenting ++// row reduction and augmentation phases require. Any `v_init` gives the ++// same optimum; duals close to the optimum leave few rows to augment. ++// ++// Kept in its own translation unit so that lap()'s code generation, which is ++// sensitive to layout, is unaffected. Final column duals are left in ++// scratch.v. ++ ++#include "lap.h" ++#include ++ ++namespace { ++inline bool nontrivially_less_than(cost a, cost b) noexcept { ++ return a + ((a > ROUND_PRECISION) ? 8 : 0) < b; ++} ++} ++ ++cost lap_warm(const lap_row dim, CostMatrix &input_cost, ++ std::vector &rowsol, std::vector &colsol, ++ LapScratch &scratch, const cost *v_init) { ++ lap_row num_free = 0; ++ scratch.ensure(dim); ++ auto &v = scratch.v; ++ auto &freeunassigned = scratch.freeunassigned; ++ const cost *__restrict__ v_ptr = v.data(); ++ ++ // WARM START ++ for (lap_col j = 0; j < dim; ++j) { ++ v[j] = v_init[j]; ++ colsol[j] = -1; ++ } ++ for (lap_row i = 0; i < dim; ++i) { ++ const cost *__restrict__ row_i = input_cost.row(i); ++ cost best = row_i[0] - v_ptr[0]; ++ lap_col best_j = 0; ++ for (lap_col j = 1; j < dim; ++j) { ++ const cost h = row_i[j] - v_ptr[j]; ++ if (h < best) { ++ best = h; ++ best_j = j; ++ } ++ } ++ if (colsol[best_j] < 0) { ++ colsol[best_j] = i; ++ rowsol[i] = best_j; ++ } else { ++ freeunassigned[num_free++] = i; ++ } ++ } ++ ++ // AUGMENTING ROW REDUCTION ++ int loopcnt = 0; ++ do { ++ ++loopcnt; ++ lap_row previous_num_free = num_free; ++ num_free = 0; ++ lap_row k = 0; ++ while (k < previous_num_free) { ++ const lap_row i = freeunassigned[k++]; ++ const auto [umin, usubmin, min_idx, j2] = ++ input_cost.findRowSubmin(&i, v); ++ lap_col j1 = min_idx; ++ lap_row i0 = colsol[j1]; ++ const bool strictly_less = nontrivially_less_than(umin, usubmin); ++ if (strictly_less) { ++ v[j1] -= (usubmin - umin); ++ } else if (i0 > -1) { ++ j1 = j2; ++ i0 = colsol[j2]; ++ } ++ rowsol[i] = j1; ++ colsol[j1] = i; ++ if (i0 > -1) { ++ if (strictly_less) { ++ freeunassigned[--k] = i0; ++ } else { ++ freeunassigned[num_free++] = i0; ++ } ++ } ++ } ++ } while (loopcnt < 2); ++ ++ // AUGMENT SOLUTION for each free row. ++ cost *__restrict__ d_ptr = scratch.d.data(); ++ lap_row *__restrict__ pred_ptr = scratch.predecessor.data(); ++ lap_col *__restrict__ cl_ptr = scratch.col_list.data(); ++ ++ for (lap_row f = 0; f < num_free; ++f) { ++ bool unassignedfound = false; ++ lap_row free_row = freeunassigned[f]; ++ const cost *__restrict__ free_row_cost = input_cost.row(free_row); ++ lap_col endofpath = 0; ++ lap_col last = 0; ++ lap_row i; ++ lap_col j1; ++ ++ for (lap_col j = 0; j < dim; ++j) { ++ d_ptr[j] = free_row_cost[j] - v_ptr[j]; ++ pred_ptr[j] = free_row; ++ cl_ptr[j] = j; ++ } ++ ++ cost min = 0; ++ lap_col low = 0; ++ lap_col up = 0; ++ ++ do { ++ if (up == low) { ++ last = low - 1; ++ min = d_ptr[cl_ptr[up++]]; ++ for (lap_dim k = up; k < dim; ++k) { ++ const lap_col j = cl_ptr[k]; ++ const cost h = d_ptr[j]; ++ if (h <= min) { ++ if (h < min) { ++ up = low; ++ min = h; ++ } ++ cl_ptr[k] = cl_ptr[up]; ++ cl_ptr[up++] = j; ++ } ++ } ++ for (lap_dim k = low; k < up; ++k) { ++ if (colsol[cl_ptr[k]] < 0) { ++ endofpath = cl_ptr[k]; ++ unassignedfound = true; ++ break; ++ } ++ } ++ } ++ ++ if (!unassignedfound) { ++ j1 = cl_ptr[low++]; ++ i = colsol[j1]; ++ const cost *__restrict__ row_i = input_cost.row(i); ++ const cost h = row_i[j1] - v_ptr[j1] - min; ++ for (lap_dim k = up; k < dim; ++k) { ++ const lap_col j = cl_ptr[k]; ++ cost v2 = row_i[j] - v_ptr[j] - h; ++ if (v2 < d_ptr[j]) { ++ pred_ptr[j] = i; ++ if (v2 == min) { ++ if (colsol[j] < 0) { ++ endofpath = j; ++ unassignedfound = true; ++ break; ++ } else { ++ cl_ptr[k] = cl_ptr[up]; ++ cl_ptr[up++] = j; ++ } ++ } ++ d_ptr[j] = v2; ++ } ++ } ++ } ++ } while (!unassignedfound); ++ ++ for (lap_dim k = 0; k <= last; ++k) { ++ j1 = cl_ptr[k]; ++ v[j1] += d_ptr[j1] - min; ++ } ++ ++ do { ++ i = pred_ptr[endofpath]; ++ colsol[endofpath] = i; ++ j1 = endofpath; ++ endofpath = rowsol[i]; ++ rowsol[i] = j1; ++ } while (i != free_row); ++ } ++ ++ cost lapcost = 0; ++ for (lap_dim i = 0; i < dim; ++i) { ++ lapcost += input_cost(i, rowsol[i]); ++ } ++ return lapcost; ++} +diff --git a/src/spr.cpp b/src/spr.cpp +index 7651639..e52d3f8 100644 +--- a/src/spr.cpp ++++ b/src/spr.cpp +@@ -5,6 +5,7 @@ + #include // for root_on_node + #include // for SplitList, count_bits + #include "reduce_tree.h" ++#include "col_prices.h" + #include // for log2 + #include // for sort, lower_bound, equal + #include +@@ -248,7 +249,10 @@ List keep_and_reduce(const List& tree1, + extern double mutual_clustering_bounded_score(const RawMatrix &x, + const RawMatrix &y, + const int32 n_tip, +- const double min_score); ++ const double min_score, ++ const ColPrices *prices); ++extern std::vector mutual_clustering_col_prices( ++ const RawMatrix &x, const RawMatrix &y, const int32 n_tips); + + namespace { + +@@ -261,9 +265,45 @@ struct TreeSplits { + double entropy; // As ClusteringEntropy() + }; + ++using price_t = int_fast64_t; ++ ++// Prices (optimal column duals) for matching the splits of a pair of trees ++// whose candidate leaves are being scored. Each candidate's reduced trees ++// are these trees restricted to fewer leaves, so each candidate split is the ++// restriction of a split here, whose price is near-optimal for it. ++struct ParentPrices { ++ int n_tip = 0; ++ int n_word = 0; // 64-bit words per split ++ int n_split = 0; ++ std::vector bits; // Splits of tree 2, n_split * n_word ++ std::vector price; ++ std::vector label; // CHARSXP of each leaf ++}; ++ ++// Position of each of `labels` among the parent's leaves; false unless ++// `labels` is a subsequence of the parent's labels. ++bool MapLeaves(const ParentPrices &pp, const CharacterVector &labels, ++ std::vector &map) { ++ const int n_tip = labels.size(); ++ map.resize(n_tip); ++ int j = 0; ++ for (int t = 0; t != n_tip; ++t) { ++ const SEXP lab = STRING_ELT(labels, t); ++ while (j < pp.n_tip && pp.label[j] != lab) ++j; ++ if (j == pp.n_tip) return false; ++ map[t] = j++; ++ } ++ return n_tip > 0; ++} ++ + // As TreeTools::as.Splits(): splits in node order, omitting the root and + // its non-tip child. Edges must be in postorder. +-TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { ++// If `leaf_index` is given, also writes each split to `indexed` as ++// `n_word` 64-bit words in which leaf t sets bit leaf_index[t]. ++TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip, ++ const std::vector *leaf_index = nullptr, ++ const int n_word = 0, ++ std::vector *indexed = nullptr) { + const int n_edge = edge.nrow(); + const int n_vertex = n_edge + 1; + const int n_bin = (n_tip - 1) / SPLIT_BIN + 1; +@@ -285,6 +325,22 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { + } + } + ++ static thread_local std::vector below_ix; ++ if (leaf_index) { ++ below_ix.assign(size_t(n_vertex) * n_word, 0); ++ for (int t = 0; t != n_tip; ++t) { ++ const int ix = (*leaf_index)[t]; ++ below_ix[size_t(t) * n_word + (ix >> 6)] = uint64_t(1) << (ix & 63); ++ } ++ for (int e = 0; e != n_edge; ++e) { ++ uint64_t *p = below_ix.data() + size_t(parent[e] - 1) * n_word; ++ const uint64_t *c = below_ix.data() + size_t(child[e] - 1) * n_word; ++ for (int w = 0; w != n_word; ++w) { ++ p[w] |= c[w]; ++ } ++ } ++ } ++ + const int root_node = parent[n_edge - 1]; + int root_child = child[n_edge - 1]; + int root_children = 1; +@@ -308,6 +364,13 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { + } + const int n_split = rows.size(); + RawMatrix ret(n_split, n_bin); ++ if (leaf_index) { ++ indexed->resize(size_t(n_split) * n_word); ++ for (int r = 0; r != n_split; ++r) { ++ std::copy_n(below_ix.data() + size_t(rows[r]) * n_word, n_word, ++ indexed->data() + size_t(r) * n_word); ++ } ++ } + // Sum of n * H(split), in bits + const double n_lg2_n = n_tip * std::log2(double(n_tip)); + double n_entropy = 0; +@@ -325,23 +388,129 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { + return {ret, n_entropy / n_tip}; + } + ++// Prices for the splits of tree 2 of a candidate, from those of its parent. ++struct CandidatePrices { ++ const ParentPrices *pp; ++ const std::vector *map; // Parent position of each candidate leaf ++ const std::vector *bits; // Candidate splits, parent leaf indexing ++ int n_split; ++ std::vector price; ++ ++ // Each candidate split takes the least price among the parent splits that ++ // restrict to it (a split of the candidate tree always has one, as the ++ // candidate tree is the parent tree restricted to fewer leaves). ++ static const price_t *Get(void *context) { ++ CandidatePrices &cp = *static_cast(context); ++ const ParentPrices &pp = *cp.pp; ++ const std::vector &map = *cp.map; ++ const int n_word = pp.n_word; ++ static thread_local std::vector mask, keys; ++ static thread_local std::vector key_price; ++ static thread_local std::vector slot; ++ mask.assign(n_word, 0); ++ for (const int ix : map) { ++ mask[ix >> 6] |= uint64_t(1) << (ix & 63); ++ } ++ // Polarize so that the candidate's first leaf lies outside each split ++ const int ref_word = map[0] >> 6; ++ const uint64_t ref_bit = uint64_t(1) << (map[0] & 63); ++ int capacity = 4; ++ while (capacity < 2 * pp.n_split) capacity <<= 1; ++ slot.assign(capacity, -1); ++ keys.resize(size_t(pp.n_split) * n_word); ++ key_price.resize(pp.n_split); ++ const auto Hash = [n_word, capacity](const uint64_t *k) { ++ uint64_t h = 0x9E3779B97F4A7C15ULL; ++ for (int w = 0; w != n_word; ++w) { ++ h ^= k[w]; ++ h *= 0xBF58476D1CE4E5B9ULL; ++ h ^= h >> 31; ++ } ++ return int(h & uint64_t(capacity - 1)); ++ }; ++ int n_key = 0; ++ for (int q = 0; q != pp.n_split; ++q) { ++ const uint64_t *src = &pp.bits[size_t(q) * n_word]; ++ uint64_t *k = &keys[size_t(n_key) * n_word]; ++ const uint64_t flip = (src[ref_word] & ref_bit) ? ~uint64_t(0) : 0; ++ for (int w = 0; w != n_word; ++w) k[w] = (src[w] ^ flip) & mask[w]; ++ for (int h = Hash(k); ; h = (h + 1) & (capacity - 1)) { ++ const int e = slot[h]; ++ if (e < 0) { ++ slot[h] = n_key; ++ key_price[n_key++] = pp.price[q]; ++ break; ++ } ++ if (std::equal(k, k + n_word, &keys[size_t(e) * n_word])) { ++ key_price[e] = std::min(key_price[e], pp.price[q]); ++ break; ++ } ++ } ++ } ++ cp.price.assign(cp.n_split, -1); ++ std::vector key(n_word); ++ int n_found = 0; ++ long double found_sum = 0; ++ for (int i = 0; i != cp.n_split; ++i) { ++ const uint64_t *src = cp.bits->data() + size_t(i) * n_word; ++ const uint64_t flip = (src[ref_word] & ref_bit) ? ~uint64_t(0) : 0; ++ for (int w = 0; w != n_word; ++w) key[w] = (src[w] ^ flip) & mask[w]; ++ for (int h = Hash(key.data()); slot[h] >= 0; ++ h = (h + 1) & (capacity - 1)) { ++ const int e = slot[h]; ++ if (std::equal(key.begin(), key.end(), &keys[size_t(e) * n_word])) { ++ cp.price[i] = key_price[e]; ++ found_sum += key_price[e]; ++ ++n_found; ++ break; ++ } ++ } ++ } ++ if (n_found == 0) return nullptr; // #nocov ++ if (n_found < cp.n_split) { ++ const price_t fill = price_t(found_sum / n_found); // #nocov ++ for (price_t &v : cp.price) if (v < 0) v = fill; // #nocov ++ } ++ return cp.price.data(); ++ } ++}; ++ + // As ClusteringInfoDist(tree1, tree2), for trees with identical tip labels + // and edges in postorder; but if the distance exceeds `cutoff` by more than + // numerical noise, may instead return a lower bound that exceeds `cutoff`. ++// `pp` (may be null) holds prices for a parent pair of trees; `map` gives the ++// parent position of each leaf, and `bits2` the splits of tree 2 in parent ++// leaf indexing. + double BoundedCID(const TreeSplits &splits1, const TreeSplits &splits2, +- const int n_tip, const double cutoff) { ++ const int n_tip, const double cutoff, ++ const ParentPrices *pp = nullptr, ++ const std::vector *map = nullptr, ++ const std::vector *bits2 = nullptr) { + const double info = splits1.entropy + splits2.entropy; + const double tol = std::max(1.0, info) * SQRT_EPS; +- const double mci = n_tip < 4 ? 0 : +- mutual_clustering_bounded_score(splits1.bits, splits2.bits, n_tip, +- (info - cutoff) / 2 - tol); ++ double mci = 0; ++ if (n_tip >= 4) { ++ CandidatePrices cp{pp, map, bits2, int(splits2.bits.nrow()), {}}; ++ const ColPrices prices{&CandidatePrices::Get, &cp}; ++ mci = mutual_clustering_bounded_score(splits1.bits, splits2.bits, n_tip, ++ (info - cutoff) / 2 - tol, ++ pp ? &prices : nullptr); ++ } + const double ret = info - mci - mci; + return ret < tol ? 0 : ret; + } + +-double BoundedCID(const List &tree1, const List &tree2, const double cutoff) { ++double BoundedCID(const List &tree1, const List &tree2, const double cutoff, ++ const ParentPrices *pp = nullptr) { + const CharacterVector labels = tree1["tip.label"]; + const int n_tip = labels.size(); ++ std::vector map; ++ std::vector bits2; ++ if (pp && MapLeaves(*pp, labels, map)) { ++ return BoundedCID(GetSplits(tree1["edge"], n_tip), ++ GetSplits(tree2["edge"], n_tip, &map, pp->n_word, &bits2), ++ n_tip, cutoff, pp, &map, &bits2); ++ } + return BoundedCID(GetSplits(tree1["edge"], n_tip), + GetSplits(tree2["edge"], n_tip), n_tip, cutoff); + } +@@ -425,7 +594,8 @@ List KeepAndReduce(const List &tree1, const List &tree2, + // `pre1` and `pre2` are the trees' edges in preorder. + double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + const CharacterVector &labels, +- const LogicalVector &keep, const double cutoff) { ++ const LogicalVector &keep, const double cutoff, ++ const ParentPrices *pp) { + const List outcome = KeepAndReduce(pre1, pre2, labels, keep); + if (Rf_isNull(outcome[0])) { + return R_NegInf; +@@ -434,12 +604,19 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + const List out2 = outcome[1]; + const CharacterVector out_labels = out1["tip.label"]; + const int n_tip = out_labels.size(); ++ std::vector map; ++ std::vector bits2; ++ if (pp && !MapLeaves(*pp, out_labels, map)) { ++ pp = nullptr; // #nocov ++ } + const TreeSplits splits1 = GetSplits(out1["edge"], n_tip); +- const TreeSplits splits2 = GetSplits(out2["edge"], n_tip); ++ const TreeSplits splits2 = pp ? ++ GetSplits(out2["edge"], n_tip, &map, pp->n_word, &bits2) : ++ GetSplits(out2["edge"], n_tip); + const int first_match = FirstMatchingSplit(splits1.bits, splits2.bits, + n_tip); + if (first_match < 0) { +- return BoundedCID(splits1, splits2, n_tip, cutoff); ++ return BoundedCID(splits1, splits2, n_tip, cutoff, pp, &map, &bits2); + } + + // Divide at the shared split, anchoring each half to the other +@@ -460,7 +637,7 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + if (Rf_isNull(sub1[0])) { + Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov + } +- const double score1 = BoundedCID(sub1[0], sub1[1], cutoff); ++ const double score1 = BoundedCID(sub1[0], sub1[1], cutoff, pp); + if (score1 > cutoff) { + return score1; + } +@@ -468,7 +645,7 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, + if (Rf_isNull(sub2[0])) { + Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov + } +- return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1); ++ return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1, pp); + } + + } +@@ -484,6 +661,27 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { + const int n_tip = labels.size(); + const IntegerMatrix pre1 = reverse(tree1["edge"]); + const IntegerMatrix pre2 = reverse(tree2["edge"]); ++ ++ // Prices from the full trees, to bound and warm-start each candidate ++ ParentPrices pp; ++ pp.n_tip = n_tip; ++ pp.n_word = (n_tip + 63) / 64; ++ pp.label.resize(n_tip); ++ for (int t = 0; t != n_tip; ++t) { ++ pp.label[t] = STRING_ELT(labels, t); ++ } ++ std::vector identity(n_tip); ++ for (int t = 0; t != n_tip; ++t) { ++ identity[t] = t; ++ } ++ { ++ const TreeSplits full1 = GetSplits(tree1["edge"], n_tip); ++ const TreeSplits full2 = GetSplits(tree2["edge"], n_tip, &identity, ++ pp.n_word, &pp.bits); ++ pp.price = mutual_clustering_col_prices(full1.bits, full2.bits, n_tip); ++ pp.n_split = pp.price.size(); ++ } ++ const ParentPrices *pp_ptr = pp.n_split ? &pp : nullptr; + NumericVector scores(n_tip); + double best = R_PosInf; + bool all_finite = true; +@@ -491,7 +689,7 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { + for (int i = 0; i != n_tip; ++i) { + Rcpp::checkUserInterrupt(); + keep[i] = false; +- scores[i] = ScoreWithout(pre1, pre2, labels, keep, best); ++ scores[i] = ScoreWithout(pre1, pre2, labels, keep, best, pp_ptr); + keep[i] = true; + if (!R_FINITE(scores[i])) { + all_finite = false; +@@ -512,7 +710,8 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { + keep[i] = false; + for (int j = i + 1; j != n_tip; ++j) { + keep[j] = false; +- pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best); ++ pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best, ++ pp_ptr); + keep[j] = true; + if (!R_FINITE(pair_scores[pair])) { + return List::create(Named("scores") = scores, +diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp +index 0dc6eae..aa0b270 100644 +--- a/src/tree_distances.cpp ++++ b/src/tree_distances.cpp +@@ -5,6 +5,7 @@ + #include /* for unique_ptr, make_unique */ + #include + #include "tree_distances.h" ++#include "col_prices.h" + + using Rcpp::_; + using Rcpp::IntegerVector; +@@ -403,10 +404,26 @@ List msi_distance(const RawMatrix &x, const RawMatrix &y, const int32 n_tips) { + + } + ++// Warm-started LAP; see lap_warm.cpp. ++extern cost lap_warm(lap_row dim, CostMatrix &input_cost, ++ std::vector &rowsol, ++ std::vector &colsol, LapScratch &scratch, ++ const cost *v_init); ++ + // With Bounded = true, gives up once the score is known to fall below + // min_score. It then returns an upper bound on the score, below min_score, + // and an empty matching. Scores at or above min_score are computed exactly as + // with Bounded = false. ++// ++// Bounded only: `prices` (may be null) can supply a price v_j for each split ++// j of `y`, typically the optimal column duals of a closely related problem. ++// When both trees have the same number of splits, the matching is a perfect ++// assignment, and for any prices ++// sum_j v_j + sum_i min_j (c_ij - v_j) <= minimum assignment cost, ++// because u_i = min_j (c_ij - v_j) makes (u, v) feasible for the dual. This ++// gives an upper bound on the score once the cost matrix is filled. If that ++// bound does not reject the candidate, the prices warm-start the LAP. ++// Prices affect only speed, never the result. + // Reusable cost matrices for mutual_clustering(). + static cost_matrix &BoundedPool(const int i, const std::size_t dim) { + static thread_local cost_matrix pool[2]; +@@ -416,7 +433,8 @@ static cost_matrix &BoundedPool(const int i, const std::size_t dim) { + + template + List mutual_clustering(const RawMatrix &x, const RawMatrix &y, +- const int32 n_tips, const double min_score = 0) { ++ const int32 n_tips, const double min_score = 0, ++ const ColPrices *prices = nullptr) { + const SplitList a(x); + const SplitList b(y); + const bool a_has_more_splits = (a.n_splits > b.n_splits); +@@ -549,10 +567,64 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + + const split_int lap_dim = most_splits - exact_matches; + ASSERT(lap_dim > 0); ++ ++ // Bounded: dual bound from prices, on the problem that remains once exact ++ // matches are removed. Each row minimum ranges over every column, so it ++ // cannot exceed the minimum over the remaining columns. Prices lie in ++ // [0, BIG] and costs in [0, BIG], so |bound| <= 2 * lap_dim * BIG, which ++ // fits in a cost when most_splits <= SL_MAX_SPLITS / 2. ++ [[maybe_unused]] const cost *col_price = nullptr; ++ [[maybe_unused]] static thread_local std::vector v_init; ++ if constexpr (Bounded) { ++ if (prices && a.n_splits == b.n_splits && ++ most_splits <= SL_MAX_SPLITS / 2) { ++ col_price = prices->get(prices->context); ++ } ++ if (col_price) { ++ cost lb = 0; ++ v_init.clear(); ++ for (split_int bi = 0; bi < b.n_splits; ++bi) { ++ if (!b_match[bi]) { ++ lb += col_price[bi]; ++ v_init.push_back(col_price[bi]); ++ } ++ } ++ for (split_int ai = 0; ai < a.n_splits; ++ai) { ++ if (a_match[ai]) continue; ++ const cost *row = score.row(ai); ++ cost row_min = row[0] - col_price[0]; ++ for (split_int bi = 1; bi < b.n_splits; ++bi) { ++ const cost reduced = row[bi] - col_price[bi]; ++ if (reduced < row_min) row_min = reduced; ++ } ++ lb += row_min; ++ } ++ const cost lap_ub = max_score * lap_dim - lb; ++ const double bound = exact_matches ? ++ static_cast(lap_ub) * over_max_score + ++ exact_match_score / n_tips : ++ static_cast(lap_ub) / max_score; ++ if (bound < min_score) { ++ return List::create(Named("score") = bound, ++ _["matching"] = IntegerVector(0)); ++ } ++ } ++ } ++ + std::vector rowsol; + std::vector colsol; + resize_uninitialized(rowsol, lap_dim); + resize_uninitialized(colsol, lap_dim); ++ [[maybe_unused]] static thread_local LapScratch warm_scratch; ++ const auto Solve = [&](cost_matrix &m) -> cost { ++ if constexpr (Bounded) { ++ if (col_price) { ++ return lap_warm(lap_dim, m, rowsol, colsol, warm_scratch, ++ v_init.data()); ++ } ++ } ++ return lap(lap_dim, m, rowsol, colsol, !Bounded); ++ }; + cost_matrix local_small(Bounded ? 0 : lap_dim); + cost_matrix &small_score = Bounded ? BoundedPool(1, lap_dim) : local_small; + +@@ -572,7 +644,7 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + small_score.padAfterRow(lap_dim - b_extra_splits, max_score); + + const double lap_score = static_cast((max_score * lap_dim) - +- lap(lap_dim, small_score, rowsol, colsol, !Bounded)) * over_max_score; ++ Solve(small_score)) * over_max_score; + const double final_score = lap_score + (exact_match_score / n_tips); + + std::unique_ptr lap_decode = std::make_unique(lap_dim); +@@ -610,7 +682,7 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + } + + const double final_score = static_cast( +- (max_score * lap_dim) - lap(lap_dim, score, rowsol, colsol, !Bounded) ++ (max_score * lap_dim) - Solve(score) + ) / max_score; + + std::vector final_matching; +@@ -770,9 +842,71 @@ List cpp_mutual_clustering(const RawMatrix &x, const RawMatrix &y, + // min_score. Used to discard candidates early in SPRDist(method = "Rogue"). + double mutual_clustering_bounded_score(const RawMatrix &x, const RawMatrix &y, + const int32 n_tip, +- const double min_score) { ++ const double min_score, ++ const ColPrices *prices) { + return Rcpp::as( +- mutual_clustering(x, y, n_tip, min_score)["score"]); ++ mutual_clustering(x, y, n_tip, min_score, prices)["score"]); ++} ++ ++// Optimal column duals of the mutual clustering LAP between the splits `x` ++// and `y` (identical splits included), one per split in `y`, in cost units, ++// shifted so that the least is zero and capped at BIG: prices for ++// mutual_clustering() on related problems. Empty if there is nothing ++// to match. ++std::vector mutual_clustering_col_prices(const RawMatrix &x, ++ const RawMatrix &y, ++ const int32 n_tips) { ++ const SplitList a(x); ++ const SplitList b(y); ++ const split_int dim = std::max(a.n_splits, b.n_splits); ++ if (dim == 0 || n_tips < 4 || b.n_splits == 0) { ++ return {}; ++ } ++ constexpr cost max_score = BIG; ++ const double max_over_tips = static_cast(max_score) / n_tips; ++ const double lg2_n = lg2_lookup(n_tips); ++ cost_matrix score(dim); ++ for (split_int ai = 0; ai < a.n_splits; ++ai) { ++ const split_int na = a.in_split[ai]; ++ const split_int nA = n_tips - na; ++ const double offset_a = lg2_n - lg2_lookup(na); ++ const double offset_A = lg2_n - lg2_lookup(nA); ++ for (split_int bi = 0; bi < b.n_splits; ++bi) { ++ split_int a_and_b = 0; ++ for (split_int bin = 0; bin < a.n_bins; ++bin) { ++ a_and_b += count_bits(a.state[ai][bin] & b.state[bi][bin]); ++ } ++ const split_int nb = b.in_split[bi]; ++ const split_int nB = n_tips - nb; ++ const split_int a_and_B = na - a_and_b; ++ const split_int A_and_b = nb - a_and_b; ++ const split_int A_and_B = nA - A_and_b; ++ if (a_and_b == A_and_b && a_and_b == a_and_B && a_and_b == A_and_B) { ++ score(ai, bi) = max_score; ++ } else { ++ const double lg2_nb = lg2_lookup(nb); ++ const double lg2_nB = lg2_lookup(nB); ++ const double ic_sum = ++ a_and_b * (lg2_lookup(a_and_b) + offset_a - lg2_nb) + ++ a_and_B * (lg2_lookup(a_and_B) + offset_a - lg2_nB) + ++ A_and_b * (lg2_lookup(A_and_b) + offset_A - lg2_nb) + ++ A_and_B * (lg2_lookup(A_and_B) + offset_A - lg2_nB); ++ score(ai, bi) = max_score - static_cast(ic_sum * max_over_tips); ++ } ++ } ++ score.padRowAfterCol(ai, b.n_splits, max_score); ++ } ++ score.padAfterRow(a.n_splits, max_score); ++ std::vector rowsol(dim); ++ std::vector colsol(dim); ++ LapScratch scratch; ++ lap(dim, score, rowsol, colsol, false, scratch); ++ std::vector ret(scratch.v.begin(), scratch.v.begin() + b.n_splits); ++ const cost least = *std::min_element(ret.begin(), ret.end()); ++ for (cost &v : ret) { ++ v = std::min(v - least, max_score); ++ } ++ return ret; + } + + // [[Rcpp::export]] +@@ -782,7 +916,7 @@ double cpp_mutual_clustering_bounded(const RawMatrix &x, const RawMatrix &y, + ASSERT(x.cols() == y.cols() && "Input splits must address same number of tips."); + const int32 n_tip = static_cast(nTip[0]); + TreeDist::check_ntip(n_tip); +- return mutual_clustering_bounded_score(x, y, n_tip, minScore); ++ return mutual_clustering_bounded_score(x, y, n_tip, minScore, nullptr); + } + + // [[Rcpp::export]] diff --git a/dev/rogue_bound/post_arr_exit.patch b/dev/rogue_bound/post_arr_exit.patch new file mode 100644 index 00000000..d6de8ed5 --- /dev/null +++ b/dev/rogue_bound/post_arr_exit.patch @@ -0,0 +1,318 @@ +diff --git a/src/lap_arr.cpp b/src/lap_arr.cpp +new file mode 100644 +index 0000000..e59e672 +--- /dev/null ++++ b/src/lap_arr.cpp +@@ -0,0 +1,240 @@ ++// Cold JV LAP (copy of TreeDist::lap()) that, after augmenting row ++// reduction, abandons the solve if the dual bound exceeds `reject_above`. ++// Scratch research only. ++#include "lap.h" ++ ++namespace { ++inline bool nontrivially_less_than(cost a, cost b) noexcept { ++ return a + ((a > ROUND_PRECISION) ? 8 : 0) < b; ++} ++} ++#define TREEDIST_CHECK_INTERRUPT() ((void)0) ++using namespace TreeDist; ++cost lap_arr_bounded(const lap_row dim, ++ CostMatrix& input_cost, ++ std::vector& rowsol, ++ std::vector& colsol, ++ LapScratch& scratch, ++ const cost reject_above, bool &rejected) ++{ ++ constexpr bool allow_interrupt = false; ++ rejected = false; ++ lap_row num_free = 0; ++ scratch.ensure(dim); ++ auto& v = scratch.v; ++ auto& matches = scratch.matches; ++ std::fill(matches.begin(), matches.begin() + dim, 0); ++ const cost* __restrict__ v_ptr = v.data(); ++ ++ // COLUMN REDUCTION ++ for (lap_col j = dim; j--; ) { ++ const auto [min, imin] = input_cost.findColMin(j); ++ v[j] = min; ++ ++matches[imin]; ++ ++ if (matches[imin] == 1) { ++ rowsol[imin] = j; ++ colsol[j] = imin; ++ } else if (v_ptr[j] < v_ptr[rowsol[imin]]) { ++ const lap_col j1 = rowsol[imin]; ++ rowsol[imin] = j; ++ colsol[j] = imin; ++ colsol[j1] = -1; ++ } else { ++ colsol[j] = -1; ++ } ++ } ++ ++ // REDUCTION TRANSFER ++ auto& freeunassigned = scratch.freeunassigned; ++ ++ for (lap_row i = 0; i < dim; ++i) { ++ if (matches[i] == 0) { ++ freeunassigned[num_free++] = i; ++ } else if (matches[i] == 1) { ++ const lap_col j1 = rowsol[i]; ++ const cost* row_i = input_cost.row(i); ++ cost min_cost; ++ if (j1 == 0) { ++ min_cost = row_i[1] - v_ptr[1]; ++ for (lap_col j = 2; j < dim; ++j) { ++ const cost reduced_cost = row_i[j] - v_ptr[j]; ++ if (reduced_cost < min_cost) { ++ min_cost = reduced_cost; ++ } ++ } ++ } else { ++ min_cost = row_i[0] - v_ptr[0]; ++ for (lap_col j = 1; j < dim; ++j) { ++ if (j == j1) continue; ++ const cost reduced_cost = row_i[j] - v_ptr[j]; ++ if (reduced_cost < min_cost) { ++ min_cost = reduced_cost; ++ } ++ } ++ } ++ v[j1] -= min_cost; ++ } ++ } ++ ++ // AUGMENTING ROW REDUCTION ++ int loopcnt = 0; ++ ++ do { ++ ++loopcnt; ++ ++ lap_row previous_num_free = num_free; ++ num_free = 0; ++ lap_row k = 0; ++ while (k < previous_num_free) { ++ const lap_row i = freeunassigned[k++]; ++ const auto [umin, usubmin, min_idx, j2] = ++ input_cost.findRowSubmin(&i, v); ++ lap_col j1 = min_idx; ++ ++ lap_row i0 = colsol[j1]; ++ const bool strictly_less = ++ nontrivially_less_than(umin, usubmin); ++ if (strictly_less) { ++ v[j1] -= (usubmin - umin); ++ } else if (i0 > -1) { ++ j1 = j2; ++ i0 = colsol[j2]; ++ } ++ ++ rowsol[i] = j1; ++ colsol[j1] = i; ++ ++ if (i0 > -1) { ++ if (strictly_less) { ++ freeunassigned[--k] = i0; ++ if (allow_interrupt) TREEDIST_CHECK_INTERRUPT(); ++ } else { ++ freeunassigned[num_free++] = i0; ++ } ++ } ++ } ++ } while (loopcnt < 2); ++ ++ // Dual bound after ARR: sum_j v_j + sum_i min_j (c_ij - v_j) is a lower ++ // bound on the assignment cost for any v. ++ if (num_free) { ++ cost lb = 0; ++ for (lap_col j = 0; j < dim; ++j) lb += v_ptr[j]; ++ for (lap_row i = 0; i < dim; ++i) { ++ const cost* row_i = input_cost.row(i); ++ cost m = row_i[0] - v_ptr[0]; ++ for (lap_col j = 1; j < dim; ++j) { ++ const cost h = row_i[j] - v_ptr[j]; ++ if (h < m) m = h; ++ } ++ lb += m; ++ } ++ if (lb > reject_above) { ++ rejected = true; ++ return lb; ++ } ++ } ++ ++ // AUGMENT SOLUTION for each free row. ++ // Restrict-qualified local pointers enable the compiler to avoid ++ // reloads after stores in the Dijkstra inner loop. ++ cost* __restrict__ d_ptr = scratch.d.data(); ++ lap_row* __restrict__ pred_ptr = scratch.predecessor.data(); ++ lap_col* __restrict__ cl_ptr = scratch.col_list.data(); ++ ++ for (lap_row f = 0; f < num_free; ++f) { ++ bool unassignedfound = false; ++ lap_row free_row = freeunassigned[f]; ++ const cost* __restrict__ free_row_cost = input_cost.row(free_row); ++ lap_col endofpath = 0; ++ lap_col last = 0; ++ lap_row i; ++ lap_col j1; ++ ++ for (lap_col j = 0; j < dim; ++j) { ++ d_ptr[j] = free_row_cost[j] - v_ptr[j]; ++ pred_ptr[j] = free_row; ++ cl_ptr[j] = j; ++ } ++ ++ cost min = 0; ++ lap_col low = 0; ++ lap_col up = 0; ++ ++ do { ++ if (up == low) { ++ last = low - 1; ++ min = d_ptr[cl_ptr[up++]]; ++ ++ for (lap_dim k = up; k < dim; ++k) { ++ const lap_col j = cl_ptr[k]; ++ const cost h = d_ptr[j]; ++ if (h <= min) { ++ if (h < min) { ++ up = low; ++ min = h; ++ } ++ cl_ptr[k] = cl_ptr[up]; ++ cl_ptr[up++] = j; ++ } ++ } ++ for (lap_dim k = low; k < up; ++k) { ++ if (colsol[cl_ptr[k]] < 0) { ++ endofpath = cl_ptr[k]; ++ unassignedfound = true; ++ break; ++ } ++ } ++ } ++ ++ if (!unassignedfound) { ++ j1 = cl_ptr[low++]; ++ i = colsol[j1]; ++ const cost* __restrict__ row_i = input_cost.row(i); ++ const cost h = row_i[j1] - v_ptr[j1] - min; ++ ++ for (lap_dim k = up; k < dim; ++k) { ++ const lap_col j = cl_ptr[k]; ++ cost v2 = row_i[j] - v_ptr[j] - h; ++ if (v2 < d_ptr[j]) { ++ pred_ptr[j] = i; ++ if (v2 == min) { ++ if (colsol[j] < 0) { ++ endofpath = j; ++ unassignedfound = true; ++ break; ++ } else { ++ cl_ptr[k] = cl_ptr[up]; ++ cl_ptr[up++] = j; ++ } ++ } ++ d_ptr[j] = v2; ++ } ++ } ++ } ++ } while (!unassignedfound); ++ ++ for (lap_dim k = 0; k <= last; ++k) { ++ j1 = cl_ptr[k]; ++ v[j1] += d_ptr[j1] - min; ++ } ++ ++ do { ++ i = pred_ptr[endofpath]; ++ colsol[endofpath] = i; ++ j1 = endofpath; ++ endofpath = rowsol[i]; ++ rowsol[i] = j1; ++ } while (i != free_row); ++ } ++ ++ // Calculate optimal cost. ++ cost lapcost = 0; ++ for (lap_dim i = 0; i < dim; ++i) { ++ lapcost += input_cost(i, rowsol[i]); ++ } ++ ++ return lapcost; ++} ++ +diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp +index 0dc6eae..d54df48 100644 +--- a/src/tree_distances.cpp ++++ b/src/tree_distances.cpp +@@ -403,6 +403,11 @@ List msi_distance(const RawMatrix &x, const RawMatrix &y, const int32 n_tips) { + + } + ++extern cost lap_arr_bounded(const lap_row dim, CostMatrix& input_cost, ++ std::vector& rowsol, ++ std::vector& colsol, LapScratch& scratch, ++ const cost reject_above, bool &rejected); ++ + // With Bounded = true, gives up once the score is known to fall below + // min_score. It then returns an upper bound on the score, below min_score, + // and an empty matching. Scores at or above min_score are computed exactly as +@@ -553,6 +558,27 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + std::vector colsol; + resize_uninitialized(rowsol, lap_dim); + resize_uninitialized(colsol, lap_dim); ++ // Bounded: reject (in cost units) once the LAP cost is known to exceed ++ // max_score * lap_dim - (min_score - exact part) * max_score. ++ [[maybe_unused]] bool rejected = false; ++ [[maybe_unused]] cost reject_above = std::numeric_limits::max(); ++ if constexpr (Bounded) { ++ const double room = (min_score - exact_match_score / n_tips) * ++ static_cast(max_score); ++ const double limit = static_cast(max_score) * lap_dim - room; ++ // Conservative by 4096 units (~1.5e-11 bits) to absorb rounding in `limit` ++ if (limit < 9.0e18 && limit > -9.0e18) { ++ reject_above = static_cast(limit) + 4096; ++ } ++ } ++ [[maybe_unused]] static thread_local LapScratch arr_scratch; ++ const auto Solve = [&](cost_matrix &m) -> cost { ++ if constexpr (Bounded) { ++ return lap_arr_bounded(lap_dim, m, rowsol, colsol, arr_scratch, ++ reject_above, rejected); ++ } ++ return lap(lap_dim, m, rowsol, colsol, !Bounded); ++ }; + cost_matrix local_small(Bounded ? 0 : lap_dim); + cost_matrix &small_score = Bounded ? BoundedPool(1, lap_dim) : local_small; + +@@ -572,8 +598,12 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + small_score.padAfterRow(lap_dim - b_extra_splits, max_score); + + const double lap_score = static_cast((max_score * lap_dim) - +- lap(lap_dim, small_score, rowsol, colsol, !Bounded)) * over_max_score; ++ Solve(small_score)) * over_max_score; + const double final_score = lap_score + (exact_match_score / n_tips); ++ if (rejected) { ++ return List::create(Named("score") = final_score, ++ _["matching"] = IntegerVector(0)); ++ } + + std::unique_ptr lap_decode = std::make_unique(lap_dim); + split_int fuzzy_match = 0; +@@ -610,8 +640,12 @@ List mutual_clustering(const RawMatrix &x, const RawMatrix &y, + } + + const double final_score = static_cast( +- (max_score * lap_dim) - lap(lap_dim, score, rowsol, colsol, !Bounded) ++ (max_score * lap_dim) - Solve(score) + ) / max_score; ++ if (rejected) { ++ return List::create(Named("score") = final_score, ++ _["matching"] = IntegerVector(0)); ++ } + + std::vector final_matching; + final_matching.reserve(a.n_splits); diff --git a/dev/rogue_bound/profile.R b/dev/rogue_bound/profile.R new file mode 100644 index 00000000..0c6393b2 --- /dev/null +++ b/dev/rogue_bound/profile.R @@ -0,0 +1,30 @@ +# Usage: Rscript profile.R [depth] [dumpfile] +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +suppressMessages(library(TreeTools)) +wl <- readRDS("/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/workloads.rds") +pairs <- wl[[args[2]]] +if (length(args) > 4) pairs <- pairs[eval(parse(text = args[5]))] +depth <- if (length(args) > 2) as.integer(args[3]) else 1L +options(sprDepth = depth) +dumpFile <- if (length(args) > 3) args[4] else "" +Run <- function() vapply(pairs, function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +invisible(Run()) # warm up +invisible(TreeDist:::instr_get(TRUE)) +if (nzchar(dumpFile)) TreeDist:::instr_dump(dumpFile) +el <- system.time(r <- Run())[["elapsed"]] +TreeDist:::instr_dump("") +x <- TreeDist:::instr_get(TRUE) +timers <- c("total", "keepReduce", "splits", "firstMatch", "boundedCID", + "mciSetup", "mciFill", "mciCopy", "mciLap", "mciTail", "halfKeepReduce", "parent") +counts <- c("calls", "entropyExit", "rowExit", "exactAll", "lap", "lapBelow", + "whole", "half1", "half2", "candidates", "extraExit") +tk <- x[seq_along(timers)] / 2.1e9 +names(tk) <- timers +ct <- x[length(timers) + seq_along(counts)] +names(ct) <- counts +cat(sprintf("== %s %s depth %d: elapsed %.3fs, sum %g\n", basename(args[1]), args[2], depth, el, sum(r))) +cat("Seconds (C++ timers):\n"); print(round(tk, 3)) +cat("Share of elapsed:\n"); print(round(tk / el, 3)) +cat("Counts:\n"); print(ct) diff --git a/dev/rogue_bound/profile2.R b/dev/rogue_bound/profile2.R new file mode 100644 index 00000000..b42dad73 --- /dev/null +++ b/dev/rogue_bound/profile2.R @@ -0,0 +1,26 @@ +# Usage: Rscript profile2.R [depth] [subset] +# For the instrumented prototype (td_new_instr). +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +suppressMessages(library(TreeTools)) +wl <- readRDS("/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/workloads.rds") +pairs <- wl[[args[2]]] +if (length(args) > 3) pairs <- pairs[eval(parse(text = args[4]))] +depth <- if (length(args) > 2) as.integer(args[3]) else 1L +options(sprDepth = depth) +Run <- function() vapply(pairs, function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +invisible(Run()) +invisible(TreeDist:::instr_get(TRUE)) +el <- system.time(r <- Run())[["elapsed"]] +x <- TreeDist:::instr_get(TRUE) +timers <- c("total", "parentPrices", "childPrices", "keepReduce", "splits", + "firstMatch", "mciFill", "mciLap", "boundedCID") +counts <- c("calls", "rowExit", "pricedAvail", "pricedExit", "lap", "lapBelow", "exactAll") +tk <- x[seq_along(timers)] / 2.1e9 +names(tk) <- timers +ct <- x[length(timers) + seq_along(counts)] +names(ct) <- counts +cat(sprintf("== %s %s depth %d: elapsed %.3fs, sum %g\n", basename(args[1]), args[2], depth, el, sum(r))) +print(rbind(seconds = round(tk, 3), share = round(tk / el, 3))) +print(ct) diff --git a/dev/rogue_bound/profile3.R b/dev/rogue_bound/profile3.R new file mode 100644 index 00000000..aecd95c9 --- /dev/null +++ b/dev/rogue_bound/profile3.R @@ -0,0 +1,26 @@ +# Usage: Rscript profile2.R [depth] [subset] +# For the instrumented prototype (td_new_instr). +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +suppressMessages(library(TreeTools)) +wl <- readRDS("/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/workloads.rds") +pairs <- wl[[args[2]]] +if (length(args) > 3) pairs <- pairs[eval(parse(text = args[4]))] +depth <- if (length(args) > 2) as.integer(args[3]) else 1L +options(sprDepth = depth) +Run <- function() vapply(pairs, function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +invisible(Run()) +invisible(TreeDist:::instr_get(TRUE)) +el <- system.time(r <- Run())[["elapsed"]] +x <- TreeDist:::instr_get(TRUE) +timers <- c("total", "parentPrices", "getPrice", "lbPass", "keepReduce", "splits", + "firstMatch", "mciFill", "mciLap", "boundedCID") +counts <- c("calls", "rowExit", "pricedAvail", "pricedExit", "lap", "lapBelow", "exactAll", "parents") +tk <- x[seq_along(timers)] / 2.1e9 +names(tk) <- timers +ct <- x[length(timers) + seq_along(counts)] +names(ct) <- counts +cat(sprintf("== %s %s depth %d: elapsed %.3fs, sum %g\n", basename(args[1]), args[2], depth, el, sum(r))) +print(rbind(seconds = round(tk, 3), share = round(tk / el, 3))) +print(ct) diff --git a/dev/rogue_bound/prog_summary.R b/dev/rogue_bound/prog_summary.R new file mode 100644 index 00000000..4abe19fb --- /dev/null +++ b/dev/rogue_bound/prog_summary.R @@ -0,0 +1,25 @@ +# Progressive (row-by-row) priced bound: how much of the fill could be +# skipped if the parent-price bound were checked while rows are filled? +dir <- "/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/csv2" +for (f in c("rand30", "rand30d2", "rand80", "rand150", "near80", "near150")) { + x <- read.csv(file.path(dir, paste0(f, ".csv"))) + lapCalls <- x$ub_row >= x$min_score + wasted <- lapCalls & x$lap_exact < x$min_score + priced <- lapCalls & !is.na(x$ub_pv_minv) & x$ub_pv_minv < x$min_score + w <- x[priced, ] + Frac <- function(p) ifelse(is.na(p) | p > 1, 1, p) + # Fill cost ~ d^2; a reject after fraction f of rows saves (1 - f) of it. + Saved <- function(col, rows) { + sum((1 - Frac(x[[col]][rows])) * x$d[rows]^2) / sum(x$d[lapCalls]^2) + } + cat(sprintf("=== %s: LAP-stage calls %d, wasted %d, priced-bound rejects %d (%.1f%% of wasted)\n", + f, sum(lapCalls), sum(wasted), sum(priced), 100 * sum(priced) / max(1, sum(wasted)))) + cat(sprintf(" rejected before any row is filled (size profile + parent prices): %.1f%% of wasted\n", + 100 * sum(lapCalls & !is.na(x$prog_nat) & x$prog_nat == 0) / max(1, sum(wasted)))) + for (k in c("prog_nat", "prog_h", "prog_oracle")) { + pk <- w[[k]] + cat(sprintf(" %-11s rows needed among priced rejects: median %.2f, q25 %.2f, q75 %.2f; fill saved (d^2-weighted, all LAP-stage calls) %.1f%%\n", + k, median(Frac(pk)), quantile(Frac(pk), 0.25), quantile(Frac(pk), 0.75), + 100 * Saved(k, which(priced)))) + } +} diff --git a/dev/rogue_bound/run_bench.sh b/dev/rogue_bound/run_bench.sh new file mode 100644 index 00000000..4872a086 --- /dev/null +++ b/dev/rogue_bound/run_bench.sh @@ -0,0 +1,18 @@ +#!/bin/bash +# Interleaved A/B/C benchmark; each process: 1 warm-up + N timed runs. +cd /tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound +out=bench_results.txt +: > $out +for round in 1 2; do + for w in rand30 near80 near150 mid80 mid150 rand80 rand150; do + reps=5 + if [ "$w" = "rand150" ]; then reps=3; fi + for l in lib_base lib_arr lib_new2; do + Rscript bench.R $l $w $reps >> $out 2>&1 + done + done + for l in lib_base lib_arr lib_new2; do + Rscript bench.R $l rand30 3 2 >> $out 2>&1 + done +done +echo DONE >> $out diff --git a/dev/rogue_bound/scores.R b/dev/rogue_bound/scores.R new file mode 100644 index 00000000..e6cd4a22 --- /dev/null +++ b/dev/rogue_bound/scores.R @@ -0,0 +1,9 @@ +# Usage: Rscript scores.R +# Calls spr_rogue_scores() on each captured input; saves the outputs. +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +inputs <- readRDS(args[2]) +out <- lapply(inputs, function(x) TreeDist:::spr_rogue_scores(x[[1]], x[[2]], x[[3]])) +saveRDS(out, args[3]) +cat(basename(args[1]), length(out), "calls\n") diff --git a/dev/rogue_bound/summarize.R b/dev/rogue_bound/summarize.R new file mode 100644 index 00000000..0718102d --- /dev/null +++ b/dev/rogue_bound/summarize.R @@ -0,0 +1,50 @@ +# Summarize analyzer CSVs: which bounds would reject calls that currently +# need a LAP only to be rejected. +args <- commandArgs(TRUE) +dir <- "/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/csv" +files <- if (length(args)) args else sub("\\.csv$", "", list.files(dir, "\\.csv$")) +bounds <- c("ub_entropy", "ub_row", "ub_cr", "ub_rc", "ub_rt", "ub_arr", + "ub_alt1", "ub_alt3", "ub_alt10", "ub_pv_minv", "ub_pv_maxv", + "ub_pv_mean", "ub_pu", "ub_pv_alt1", "ub_pv_alt3", "ub_pv_warmarr", + "ub_prev", "ub_size") +for (f in files) { + x <- read.csv(file.path(dir, paste0(f, ".csv"))) + x$ub_size[is.na(x$ub_size)] <- Inf + for (b in bounds) x[[b]][is.na(x[[b]])] <- Inf + # Sanity: every bound must be >= the exact score + for (b in bounds) { + bad <- sum(x[[b]] < x$lap_exact - 1e-9) + if (bad) cat("!! INVALID", b, bad, "\n") + } + lapCalls <- x$ub_row >= x$min_score # not caught by the current row bound + wasted <- lapCalls & x$lap_exact < x$min_score + cat(sprintf("\n=== %s: %d calls; reach LAP %d; LAP only rejects %d (%.1f%%)\n", + f, nrow(x), sum(lapCalls), sum(wasted), 100 * mean(wasted[lapCalls]))) + w <- x[wasted, ] + rel <- (w$min_score - w$lap_exact) / w$lap_exact + cat(sprintf(" margin (min_score - MCI) for wasted LAPs: median %.4f bits (%.3f%% of MCI); q10 %.4f, q90 %.4f\n", + median(w$min_score - w$lap_exact), 100 * median(rel), + quantile(w$min_score - w$lap_exact, 0.1), quantile(w$min_score - w$lap_exact, 0.9))) + res <- t(vapply(bounds, function(b) { + c(rejects = mean(w[[b]] < w$min_score), + medExcess = median(w[[b]] - w$lap_exact), + medExcessPct = 100 * median((w[[b]] - w$lap_exact) / w$lap_exact)) + }, double(3))) + print(round(res, 4)) + if (nrow(w)) { + cat(sprintf(" augmentation phase: median free rows after ARR %.0f; bound rejects within aug %.1f%% of wasted (median frac of augs needed %.2f)\n", + median(w$aug_total), 100 * mean(w$aug_needed >= 0 & w$aug_needed < w$aug_total), + median(ifelse(w$aug_needed >= 0, w$aug_needed / pmax(1, w$aug_total), 1)))) + cat(sprintf(" frac of cols mapped from parent: median %.2f; multi-preimage cols median %.0f\n", + median(x$frac_mapped, na.rm = TRUE), median(x$n_multi, na.rm = TRUE))) + lx <- x[lapCalls, ] + cat(sprintf(" LAP work (all LAP calls): cold ticks %.0f, warm ticks %.0f (ratio %.2f); warm optimal %.4f\n", + mean(lx$cold_ticks), mean(lx$warm_ticks, na.rm = TRUE), + sum(lx$warm_ticks, na.rm = TRUE) / sum(lx$cold_ticks[!is.na(lx$warm_ticks)]), + mean(lx$warm_ok[!is.na(lx$warm_ticks) & lx$warm_ticks > 0]))) + cat(sprintf(" cold: free rows after init %.1f, after ARR %.1f, aug ops %.0f; warm: %.1f, %.1f, %.0f\n", + mean(lx$cold_free_init), mean(lx$cold_free_arr), mean(lx$cold_aug), + mean(lx$warm_free_init), mean(lx$warm_free_arr), mean(lx$warm_aug))) + cat(sprintf(" mean d %.1f\n", mean(lx$d))) + } +} diff --git a/dev/rogue_bound/verify.R b/dev/rogue_bound/verify.R new file mode 100644 index 00000000..9483878c --- /dev/null +++ b/dev/rogue_bound/verify.R @@ -0,0 +1,32 @@ +# Usage: Rscript verify.R [capture.rds] +# Runs SPRDist(method = "Rogue") on the verify pairs at depths 1 and 2. +# With a third argument, also records every spr_rogue_scores() input +# (from this lib's run) to that file. +args <- commandArgs(TRUE) +.libPaths(c("/root/Rlib", .libPaths())) +library(TreeDist, lib.loc = args[1]) +suppressMessages(library(TreeTools)) +pairs <- readRDS("/tmp/claude-0/-home-claude/37427f5a-6383-50c4-a6e2-4664447d3c94/scratchpad/bound/verify.rds") +capture <- length(args) > 2 +captured <- list() +if (capture) { + ns <- asNamespace("TreeDist") + orig <- get("spr_rogue_scores", ns) + env <- environment() + wrapped <- function(tree1, tree2, depth) { + assign("captured", c(get("captured", env), list(list(tree1, tree2, depth))), env) + orig(tree1, tree2, depth) + } + unlockBinding("spr_rogue_scores", ns) + assign("spr_rogue_scores", wrapped, ns) +} +nTip <- vapply(pairs, function(p) NTip(p[[1]]), 1) +options(sprDepth = 1) +d1 <- vapply(pairs, function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +options(sprDepth = 2) +small <- which(nTip <= 60) +d2 <- vapply(pairs[small], function(p) SPRDist(p[[1]], p[[2]], method = "Rogue"), 1) +saveRDS(list(d1 = d1, d2 = d2), args[2]) +if (capture) saveRDS(captured, args[3]) +cat(basename(args[1]), "pairs:", length(d1), "depth-2 pairs:", length(d2), + "sum d1", sum(d1), "sum d2", sum(d2), "captured", length(captured), "\n") diff --git a/dev/rogue_bound/verify.rds b/dev/rogue_bound/verify.rds new file mode 100644 index 00000000..7b106262 Binary files /dev/null and b/dev/rogue_bound/verify.rds differ diff --git a/dev/rogue_bound/workloads.rds b/dev/rogue_bound/workloads.rds new file mode 100644 index 00000000..dd0592aa Binary files /dev/null and b/dev/rogue_bound/workloads.rds differ diff --git a/dev/rogue_mci_bound_findings.md b/dev/rogue_mci_bound_findings.md new file mode 100644 index 00000000..ef5450be --- /dev/null +++ b/dev/rogue_mci_bound_findings.md @@ -0,0 +1,74 @@ +# Tighter MCI upper bounds for `SPRDist(method = "Rogue")` + +Research on TreeDist at e07f2e9 (head of #211, since merged). Scratch only, nothing pushed. The patches are in `patches/`, next to this file. + +## Recommendation + +Use **parent-dual prices**, both as a bound and as a warm start (`patches/parent_prices_warm_start.patch`, +550/−20 lines, applies cleanly to e07f2e9). + +1. At the start of each `spr_rogue_scores()` call, solve the MCI assignment once for the unreduced pair. Keep its optimal column duals v, shifted to [0, BIG]. +2. Give each candidate's reduced tree-2 split the dual of the parent split it restricts to. Where splits merge, take the least dual. +3. After the cost matrix is filled, and only if the row bound has not already exited, compute `LB = Σ_j v_j + Σ_i min_j (c_ij − v_j)` in one exact int64 pass. If the implied MCI bound is below `min_score`, return it. +4. Otherwise run JV warm-started from v (`lap_warm()`). + +**Lemma.** For any square cost matrix c and any vector v, LB(v) ≤ OPT. + +*Proof.* For any permutation σ, `Σ_i c_iσ(i) = Σ_i (c_iσ(i) − v_σ(i)) + Σ_j v_j ≥ Σ_i min_j (c_ij − v_j) + Σ_j v_j`. ∎ + +So the bound is valid whatever the prices are: poor prices only make it looser. Validity therefore does not depend on the leaf mapping or on the warm start. + +## Results + +**Share of currently wasted LAPs rejected** (LAPs that run only to end below `min_score`): + +| Bound | rand30 | rand80 | rand150 | near80 | near150 | Median excess at 80 tips (bits) | +|---|---|---|---|---|---|---| +| Entropy / size profile | 0% | 0% | 0% | 0% | 0% | 18–20 | +| Column reduction | 25% | 12% | 12% | 60% | 66% | 0.39 | +| Cold JV duals after ARR (the #174 idea) | 63% | 48% | 39% | 78% | 75% | 0.18 | +| Previous candidate's duals | 58% | 64% | 67% | 67% | 71% | 0.11 | +| **Parent column duals** | **82%** | **87%** | **88%** | **73%** | **86%** | **0.045** | + +For scale, the median gap `min_score − MCI` on these LAPs is 0.14–0.20 bits. + +**Wall-clock** (seconds, minimum of interleaved runs; output checksums identical across builds): + +| Workload | Baseline | Cold post-ARR exit | Parent prices + warm start | +|---|---|---|---| +| 30 tips, random | 0.423 | 0.380 (−10%) | 0.343 (−19%) | +| 30 tips, depth 2 | 3.659 | 3.293 (−10%) | 3.131 (−14%) | +| 80 tips, random | 3.997 | 3.442 (−14%) | 1.944 (−51%) | +| 150 tips, random | 18.209 | 15.795 (−13%) | 6.396 (−65%) | +| 80 tips, 15 moves | 0.242 | 0.222 (−8%) | 0.219 (−10%) | +| 150 tips, 30 moves | 1.307 | 1.300 (−1%) | 1.001 (−23%) | +| 80 tips, near-identical | 0.048 | 0.045 | 0.047 (noise) | +| 150 tips, near-identical | 0.053 | 0.049 | 0.055 (noise) | + +The LAPs that remain run 2.0×, 3.4× and 4.1× faster at 30, 80 and 150 tips. The LAP's share of runtime falls from 52–66% to 4–5%. The new overheads are the parent LAP (2–3%), the price mapping (2–3%) and the LB pass (1–4%). + +**Correctness:** +- SPRDist gives identical results on 345 pair runs: 210 pairs of 8–160 leaves at depth 1, and 135 at depth 2. +- All 650,676 candidate scores (from 5,391 calls) were re-scored by both builds. The minima, argmin sets and near-minimum scores were identical, and the largest difference was 7e-15. +- Warm and cold LAPs differ by at most 1 cost unit (3.5e-15 bits), which is within the solver's existing ARR tolerance. + +## Ideas that did not work + +- **Size-profile and within-tree compatibility bounds** are about 100× too loose. They cannot see the overlaps between the two trees, which are what the cells measure. +- **Tightening the bound as rows are filled** would skip at most 1.5% of the fill, because leaving even one or two rows unfilled adds more slack than the margin. +- **Coordinate ascent** from column reduction stalls after one round. +- **Extra refinement on top of parent prices** (another coordinate round, or a warm ARR check) rejects 5–10 more points, but each costs another O(d²) pass, so the net effect is within ±1%. + +## The earlier #174 verdict + +#174 concluded that a post-ARR dual exit gave no measurable gain. That was measured when MCI was under 0.2% of runtime, before scoring moved to C++. Re-tested now, a cold post-ARR exit is exact and saves 10–14% on random pairs. Parent prices beat it everywhere. + +## Next steps + +1. **Productionize the patch.** + - Add tests for the priced exit, the warm path, the mapping-failure branch, unequal split counts and the d > SL_MAX_SPLITS/2 guard. + - Add a NEWS.md entry and a version bump. + - Compute parent prices lazily, so near-identical pairs do not pay for the parent LAP. +2. **Speed up the cost-matrix fill**, which now takes 37–55% of runtime. + - 97.6% of candidates drop a single leaf. A child cell's quadrant counts are the parent cell's counts minus that leaf's quadrant, so each cell costs O(1). + - This would also skip KeepAndReduce and GetSplits (about 30% of runtime) for those candidates. + - It is the largest remaining lever. diff --git a/src/col_prices.h b/src/col_prices.h new file mode 100644 index 00000000..65b4eed7 --- /dev/null +++ b/src/col_prices.h @@ -0,0 +1,15 @@ +#ifndef TREEDIST_COL_PRICES_H_ +#define TREEDIST_COL_PRICES_H_ + +#include + +// Supplies, on demand, prices for the splits of the second tree in a bounded +// mutual clustering calculation: one per split, in LAP cost units, each in +// [0, BIG]. `get` may return null if no prices are available. Prices +// affect only speed, never results. +struct ColPrices { + const int_fast64_t *(*get)(void *context); + void *context; +}; + +#endif diff --git a/src/lap_warm.cpp b/src/lap_warm.cpp new file mode 100644 index 00000000..2b61b0b5 --- /dev/null +++ b/src/lap_warm.cpp @@ -0,0 +1,182 @@ +// Jonker-Volgenant LAP warm-started from given column duals. +// +// As TreeDist::lap() (inst/include/TreeDist/lap_impl.h), except that column +// reduction and reduction transfer are replaced by a start from the column +// duals `v_init`: each row takes its cheapest reduced-cost column if that +// column is still free, so every assigned edge is tight, as the augmenting +// row reduction and augmentation phases require. Any `v_init` gives the +// same optimum; duals close to the optimum leave few rows to augment. +// +// Kept in its own translation unit so that lap()'s code generation, which is +// sensitive to layout, is unaffected. Final column duals are left in +// scratch.v. + +#include "lap.h" +#include + +namespace { +inline bool nontrivially_less_than(cost a, cost b) noexcept { + return a + ((a > ROUND_PRECISION) ? 8 : 0) < b; +} +} + +cost lap_warm(const lap_row dim, CostMatrix &input_cost, + std::vector &rowsol, std::vector &colsol, + LapScratch &scratch, const cost *v_init) { + lap_row num_free = 0; + scratch.ensure(dim); + auto &v = scratch.v; + auto &freeunassigned = scratch.freeunassigned; + const cost *__restrict__ v_ptr = v.data(); + + // WARM START + for (lap_col j = 0; j < dim; ++j) { + v[j] = v_init[j]; + colsol[j] = -1; + } + for (lap_row i = 0; i < dim; ++i) { + const cost *__restrict__ row_i = input_cost.row(i); + cost best = row_i[0] - v_ptr[0]; + lap_col best_j = 0; + for (lap_col j = 1; j < dim; ++j) { + const cost h = row_i[j] - v_ptr[j]; + if (h < best) { + best = h; + best_j = j; + } + } + if (colsol[best_j] < 0) { + colsol[best_j] = i; + rowsol[i] = best_j; + } else { + freeunassigned[num_free++] = i; + } + } + + // AUGMENTING ROW REDUCTION + int loopcnt = 0; + do { + ++loopcnt; + lap_row previous_num_free = num_free; + num_free = 0; + lap_row k = 0; + while (k < previous_num_free) { + const lap_row i = freeunassigned[k++]; + const auto [umin, usubmin, min_idx, j2] = + input_cost.findRowSubmin(&i, v); + lap_col j1 = min_idx; + lap_row i0 = colsol[j1]; + const bool strictly_less = nontrivially_less_than(umin, usubmin); + if (strictly_less) { + v[j1] -= (usubmin - umin); + } else if (i0 > -1) { + j1 = j2; + i0 = colsol[j2]; + } + rowsol[i] = j1; + colsol[j1] = i; + if (i0 > -1) { + if (strictly_less) { + freeunassigned[--k] = i0; + } else { + freeunassigned[num_free++] = i0; + } + } + } + } while (loopcnt < 2); + + // AUGMENT SOLUTION for each free row. + cost *__restrict__ d_ptr = scratch.d.data(); + lap_row *__restrict__ pred_ptr = scratch.predecessor.data(); + lap_col *__restrict__ cl_ptr = scratch.col_list.data(); + + for (lap_row f = 0; f < num_free; ++f) { + bool unassignedfound = false; + lap_row free_row = freeunassigned[f]; + const cost *__restrict__ free_row_cost = input_cost.row(free_row); + lap_col endofpath = 0; + lap_col last = 0; + lap_row i; + lap_col j1; + + for (lap_col j = 0; j < dim; ++j) { + d_ptr[j] = free_row_cost[j] - v_ptr[j]; + pred_ptr[j] = free_row; + cl_ptr[j] = j; + } + + cost min = 0; + lap_col low = 0; + lap_col up = 0; + + do { + if (up == low) { + last = low - 1; + min = d_ptr[cl_ptr[up++]]; + for (lap_dim k = up; k < dim; ++k) { + const lap_col j = cl_ptr[k]; + const cost h = d_ptr[j]; + if (h <= min) { + if (h < min) { + up = low; + min = h; + } + cl_ptr[k] = cl_ptr[up]; + cl_ptr[up++] = j; + } + } + for (lap_dim k = low; k < up; ++k) { + if (colsol[cl_ptr[k]] < 0) { + endofpath = cl_ptr[k]; + unassignedfound = true; + break; + } + } + } + + if (!unassignedfound) { + j1 = cl_ptr[low++]; + i = colsol[j1]; + const cost *__restrict__ row_i = input_cost.row(i); + const cost h = row_i[j1] - v_ptr[j1] - min; + for (lap_dim k = up; k < dim; ++k) { + const lap_col j = cl_ptr[k]; + cost v2 = row_i[j] - v_ptr[j] - h; + if (v2 < d_ptr[j]) { + pred_ptr[j] = i; + if (v2 == min) { + if (colsol[j] < 0) { + endofpath = j; + unassignedfound = true; + break; + } else { + cl_ptr[k] = cl_ptr[up]; + cl_ptr[up++] = j; + } + } + d_ptr[j] = v2; + } + } + } + } while (!unassignedfound); + + for (lap_dim k = 0; k <= last; ++k) { + j1 = cl_ptr[k]; + v[j1] += d_ptr[j1] - min; + } + + do { + i = pred_ptr[endofpath]; + colsol[endofpath] = i; + j1 = endofpath; + endofpath = rowsol[i]; + rowsol[i] = j1; + } while (i != free_row); + } + + cost lapcost = 0; + for (lap_dim i = 0; i < dim; ++i) { + lapcost += input_cost(i, rowsol[i]); + } + return lapcost; +} diff --git a/src/spr.cpp b/src/spr.cpp index b253cd6f..387164f7 100644 --- a/src/spr.cpp +++ b/src/spr.cpp @@ -6,6 +6,7 @@ #include // for SplitList #include #include "reduce_tree.h" +#include "col_prices.h" #include // for log2 #include // for sort, lower_bound, equal #include @@ -263,7 +264,10 @@ List keep_and_reduce(const List& tree1, extern double mutual_clustering_bounded_score(const RawMatrix &x, const RawMatrix &y, const int32 n_tip, - const double min_score); + const double min_score, + const ColPrices *prices); +extern std::vector mutual_clustering_col_prices( + const RawMatrix &x, const RawMatrix &y, const int32 n_tips); namespace { @@ -276,9 +280,59 @@ struct TreeSplits { double entropy; // As ClusteringEntropy() }; +using price_t = int_fast64_t; + +// Prices (optimal column duals) for matching the splits of a pair of trees +// whose candidate leaves are being scored. Each candidate's reduced trees +// are these trees restricted to fewer leaves, so each candidate split is the +// restriction of a split here, whose price is near-optimal for it. +struct ParentPrices { + int n_tip = 0; + int n_word = 0; // 64-bit words per split + int n_split = 0; + std::vector bits; // Splits of tree 2, n_split * n_word + std::vector label; // CHARSXP of each leaf + RawMatrix splits1, splits2; // As as.Splits(), for computing prices + + // Prices are computed on first use: a parent LAP is wasted if the row + // bound rejects every candidate. Empty if there is nothing to match. + const std::vector &Price() const { + if (!ready) { + price = mutual_clustering_col_prices(splits1, splits2, n_tip); + ready = true; + } + return price; + } + +private: + mutable std::vector price; + mutable bool ready = false; +}; + +// Position of each of `labels` among the parent's leaves; false unless +// `labels` is a subsequence of the parent's labels. +bool MapLeaves(const ParentPrices &pp, const CharacterVector &labels, + std::vector &map) { + const int n_tip = labels.size(); + map.resize(n_tip); + int j = 0; + for (int t = 0; t != n_tip; ++t) { + const SEXP lab = STRING_ELT(labels, t); + while (j < pp.n_tip && pp.label[j] != lab) ++j; + if (j == pp.n_tip) return false; + map[t] = j++; + } + return n_tip > 0; +} + // As TreeTools::as.Splits(): splits in node order, omitting the root and // its non-tip child. Edges must be in postorder. -TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { +// If `leaf_index` is given, also writes each split to `indexed` as +// `n_word` 64-bit words in which leaf t sets bit leaf_index[t]. +TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip, + const std::vector *leaf_index = nullptr, + const int n_word = 0, + std::vector *indexed = nullptr) { const int n_edge = edge.nrow(); const int n_vertex = n_edge + 1; const int n_bin = (n_tip - 1) / SPLIT_BIN + 1; @@ -300,6 +354,22 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { } } + static thread_local std::vector below_ix; + if (leaf_index) { + below_ix.assign(size_t(n_vertex) * n_word, 0); + for (int t = 0; t != n_tip; ++t) { + const int ix = (*leaf_index)[t]; + below_ix[size_t(t) * n_word + (ix >> 6)] = uint64_t(1) << (ix & 63); + } + for (int e = 0; e != n_edge; ++e) { + uint64_t *p = below_ix.data() + size_t(parent[e] - 1) * n_word; + const uint64_t *c = below_ix.data() + size_t(child[e] - 1) * n_word; + for (int w = 0; w != n_word; ++w) { + p[w] |= c[w]; + } + } + } + const int root_node = parent[n_edge - 1]; int root_child = child[n_edge - 1]; int root_children = 1; @@ -323,6 +393,13 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { } const int n_split = rows.size(); RawMatrix ret(n_split, n_bin); + if (leaf_index) { + indexed->resize(size_t(n_split) * n_word); + for (int r = 0; r != n_split; ++r) { + std::copy_n(below_ix.data() + size_t(rows[r]) * n_word, n_word, + indexed->data() + size_t(r) * n_word); + } + } // Sum of n * H(split), in bits const double n_lg2_n = n_tip * std::log2(double(n_tip)); double n_entropy = 0; @@ -340,23 +417,131 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip) { return {ret, n_entropy / n_tip}; } +// Prices for the splits of tree 2 of a candidate, from those of its parent. +struct CandidatePrices { + const ParentPrices *pp; + const std::vector *map; // Parent position of each candidate leaf + const std::vector *bits; // Candidate splits, parent leaf indexing + int n_split; + std::vector price; + + // Each candidate split takes the least price among the parent splits that + // restrict to it (a split of the candidate tree always has one, as the + // candidate tree is the parent tree restricted to fewer leaves). + static const price_t *Get(void *context) { + CandidatePrices &cp = *static_cast(context); + const ParentPrices &pp = *cp.pp; + const std::vector &parent_price = pp.Price(); + if (parent_price.empty()) return nullptr; // #nocov + const std::vector &map = *cp.map; + const int n_word = pp.n_word; + static thread_local std::vector mask, keys; + static thread_local std::vector key_price; + static thread_local std::vector slot; + mask.assign(n_word, 0); + for (const int ix : map) { + mask[ix >> 6] |= uint64_t(1) << (ix & 63); + } + // Polarize so that the candidate's first leaf lies outside each split + const int ref_word = map[0] >> 6; + const uint64_t ref_bit = uint64_t(1) << (map[0] & 63); + int capacity = 4; + while (capacity < 2 * pp.n_split) capacity <<= 1; + slot.assign(capacity, -1); + keys.resize(size_t(pp.n_split) * n_word); + key_price.resize(pp.n_split); + const auto Hash = [n_word, capacity](const uint64_t *k) { + uint64_t h = 0x9E3779B97F4A7C15ULL; + for (int w = 0; w != n_word; ++w) { + h ^= k[w]; + h *= 0xBF58476D1CE4E5B9ULL; + h ^= h >> 31; + } + return int(h & uint64_t(capacity - 1)); + }; + int n_key = 0; + for (int q = 0; q != pp.n_split; ++q) { + const uint64_t *src = &pp.bits[size_t(q) * n_word]; + uint64_t *k = &keys[size_t(n_key) * n_word]; + const uint64_t flip = (src[ref_word] & ref_bit) ? ~uint64_t(0) : 0; + for (int w = 0; w != n_word; ++w) k[w] = (src[w] ^ flip) & mask[w]; + for (int h = Hash(k); ; h = (h + 1) & (capacity - 1)) { + const int e = slot[h]; + if (e < 0) { + slot[h] = n_key; + key_price[n_key++] = parent_price[q]; + break; + } + if (std::equal(k, k + n_word, &keys[size_t(e) * n_word])) { + key_price[e] = std::min(key_price[e], parent_price[q]); + break; + } + } + } + cp.price.assign(cp.n_split, -1); + std::vector key(n_word); + int n_found = 0; + long double found_sum = 0; + for (int i = 0; i != cp.n_split; ++i) { + const uint64_t *src = cp.bits->data() + size_t(i) * n_word; + const uint64_t flip = (src[ref_word] & ref_bit) ? ~uint64_t(0) : 0; + for (int w = 0; w != n_word; ++w) key[w] = (src[w] ^ flip) & mask[w]; + for (int h = Hash(key.data()); slot[h] >= 0; + h = (h + 1) & (capacity - 1)) { + const int e = slot[h]; + if (std::equal(key.begin(), key.end(), &keys[size_t(e) * n_word])) { + cp.price[i] = key_price[e]; + found_sum += key_price[e]; + ++n_found; + break; + } + } + } + if (n_found == 0) return nullptr; // #nocov + if (n_found < cp.n_split) { + const price_t fill = price_t(found_sum / n_found); // #nocov + for (price_t &v : cp.price) if (v < 0) v = fill; // #nocov + } + return cp.price.data(); + } +}; + // As ClusteringInfoDist(tree1, tree2), for trees with identical tip labels // and edges in postorder; but if the distance exceeds `cutoff` by more than // numerical noise, may instead return a lower bound that exceeds `cutoff`. +// `pp` (may be null) holds prices for a parent pair of trees; `map` gives the +// parent position of each leaf, and `bits2` the splits of tree 2 in parent +// leaf indexing. double BoundedCID(const TreeSplits &splits1, const TreeSplits &splits2, - const int n_tip, const double cutoff) { + const int n_tip, const double cutoff, + const ParentPrices *pp = nullptr, + const std::vector *map = nullptr, + const std::vector *bits2 = nullptr) { const double info = splits1.entropy + splits2.entropy; const double tol = std::max(1.0, info) * SQRT_EPS; - const double mci = n_tip < 4 ? 0 : - mutual_clustering_bounded_score(splits1.bits, splits2.bits, n_tip, - (info - cutoff) / 2 - tol); + double mci = 0; + if (n_tip >= 4) { + CandidatePrices cp{pp, map, bits2, int(splits2.bits.nrow()), {}}; + const ColPrices prices{&CandidatePrices::Get, &cp}; + mci = mutual_clustering_bounded_score(splits1.bits, splits2.bits, n_tip, + (info - cutoff) / 2 - tol, + pp ? &prices : nullptr); + } const double ret = info - mci - mci; return ret < tol ? 0 : ret; } -double BoundedCID(const List &tree1, const List &tree2, const double cutoff) { +double BoundedCID(const List &tree1, const List &tree2, const double cutoff, + const ParentPrices *pp = nullptr) { const CharacterVector labels = tree1["tip.label"]; const int n_tip = labels.size(); + std::vector map; + std::vector bits2; + if (pp && MapLeaves(*pp, labels, map)) { + return BoundedCID(GetSplits(tree1["edge"], n_tip), + GetSplits(tree2["edge"], n_tip, &map, pp->n_word, &bits2), + n_tip, cutoff, pp, &map, &bits2); + } return BoundedCID(GetSplits(tree1["edge"], n_tip), GetSplits(tree2["edge"], n_tip), n_tip, cutoff); } @@ -440,7 +625,8 @@ List KeepAndReduce(const List &tree1, const List &tree2, // `pre1` and `pre2` are the trees' edges in preorder. double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, const CharacterVector &labels, - const LogicalVector &keep, const double cutoff) { + const LogicalVector &keep, const double cutoff, + const ParentPrices *pp) { const List outcome = KeepAndReduce(pre1, pre2, labels, keep); if (Rf_isNull(outcome[0])) { return R_NegInf; @@ -449,12 +635,19 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, const List out2 = outcome[1]; const CharacterVector out_labels = out1["tip.label"]; const int n_tip = out_labels.size(); + std::vector map; + std::vector bits2; + if (pp && !MapLeaves(*pp, out_labels, map)) { + pp = nullptr; // #nocov + } const TreeSplits splits1 = GetSplits(out1["edge"], n_tip); - const TreeSplits splits2 = GetSplits(out2["edge"], n_tip); + const TreeSplits splits2 = pp ? + GetSplits(out2["edge"], n_tip, &map, pp->n_word, &bits2) : + GetSplits(out2["edge"], n_tip); const int first_match = FirstMatchingSplit(splits1.bits, splits2.bits, n_tip); if (first_match < 0) { - return BoundedCID(splits1, splits2, n_tip, cutoff); + return BoundedCID(splits1, splits2, n_tip, cutoff, pp, &map, &bits2); } // Divide at the shared split, anchoring each half to the other @@ -475,7 +668,7 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, if (Rf_isNull(sub1[0])) { Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov } - const double score1 = BoundedCID(sub1[0], sub1[1], cutoff); + const double score1 = BoundedCID(sub1[0], sub1[1], cutoff, pp); if (score1 > cutoff) { return score1; } @@ -483,7 +676,7 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, if (Rf_isNull(sub2[0])) { Rcpp::stop("Half-tree reduced to nothing; please report this bug."); // #nocov } - return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1); + return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1, pp); } } @@ -499,6 +692,28 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { const int n_tip = labels.size(); const IntegerMatrix pre1 = reverse(tree1["edge"]); const IntegerMatrix pre2 = reverse(tree2["edge"]); + + // Prices from the full trees, to bound and warm-start each candidate + ParentPrices pp; + pp.n_tip = n_tip; + pp.n_word = (n_tip + 63) / 64; + pp.label.resize(n_tip); + for (int t = 0; t != n_tip; ++t) { + pp.label[t] = STRING_ELT(labels, t); + } + std::vector identity(n_tip); + for (int t = 0; t != n_tip; ++t) { + identity[t] = t; + } + { + const TreeSplits full1 = GetSplits(tree1["edge"], n_tip); + const TreeSplits full2 = GetSplits(tree2["edge"], n_tip, &identity, + pp.n_word, &pp.bits); + pp.splits1 = full1.bits; + pp.splits2 = full2.bits; + pp.n_split = full2.bits.nrow(); + } + const ParentPrices *pp_ptr = pp.n_split ? &pp : nullptr; NumericVector scores(n_tip); double best = R_PosInf; bool all_finite = true; @@ -506,7 +721,7 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { for (int i = 0; i != n_tip; ++i) { Rcpp::checkUserInterrupt(); keep[i] = false; - scores[i] = ScoreWithout(pre1, pre2, labels, keep, best); + scores[i] = ScoreWithout(pre1, pre2, labels, keep, best, pp_ptr); keep[i] = true; if (!R_FINITE(scores[i])) { all_finite = false; @@ -527,7 +742,8 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { keep[i] = false; for (int j = i + 1; j != n_tip; ++j) { keep[j] = false; - pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best); + pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best, + pp_ptr); keep[j] = true; if (!R_FINITE(pair_scores[pair])) { return List::create(Named("scores") = scores, diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp index c1ceed17..7189bf9d 100644 --- a/src/tree_distances.cpp +++ b/src/tree_distances.cpp @@ -6,6 +6,7 @@ #include /* for unique_ptr, make_unique */ #include #include "tree_distances.h" +#include "col_prices.h" using Rcpp::_; using Rcpp::IntegerVector; @@ -435,10 +436,26 @@ List msi_distance(const RawMatrix &x, const RawMatrix &y, const int32 n_tips) { }); } +// Warm-started LAP; see lap_warm.cpp. +extern cost lap_warm(lap_row dim, CostMatrix &input_cost, + std::vector &rowsol, + std::vector &colsol, LapScratch &scratch, + const cost *v_init); + // With Bounded = true, gives up once the score is known to fall below // min_score. It then returns an upper bound on the score, below min_score, // and an empty matching. Scores at or above min_score are computed exactly as // with Bounded = false. +// +// Bounded only: `prices` (may be null) can supply a price v_j for each split +// j of `y`, typically the optimal column duals of a closely related problem. +// When both trees have the same number of splits, the matching is a perfect +// assignment, and for any prices +// sum_j v_j + sum_i min_j (c_ij - v_j) <= minimum assignment cost, +// because u_i = min_j (c_ij - v_j) makes (u, v) feasible for the dual. This +// gives an upper bound on the score once the cost matrix is filled. If that +// bound does not reject the candidate, the prices warm-start the LAP. +// Prices affect only speed, never the result. // Reusable cost matrices for mutual_clustering(). static cost_matrix &BoundedPool(const int i, const std::size_t dim) { static thread_local cost_matrix pool[2]; @@ -448,7 +465,8 @@ static cost_matrix &BoundedPool(const int i, const std::size_t dim) { template List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, - const int32 n_tips, const double min_score) { + const int32 n_tips, const double min_score, + const ColPrices *prices) { const SplitList a(x); const SplitList b(y); const bool a_has_more_splits = (a.n_splits > b.n_splits); @@ -581,10 +599,64 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, const split_int lap_dim = most_splits - exact_matches; ASSERT(lap_dim > 0); + + // Bounded: dual bound from prices, on the problem that remains once exact + // matches are removed. Each row minimum ranges over every column, so it + // cannot exceed the minimum over the remaining columns. Prices lie in + // [0, BIG] and costs in [0, BIG], so |bound| <= 2 * lap_dim * BIG, which + // fits in a cost when most_splits <= SL_MAX_SPLITS / 2. + [[maybe_unused]] const cost *col_price = nullptr; + [[maybe_unused]] static thread_local std::vector v_init; + if constexpr (Bounded) { + if (prices && a.n_splits == b.n_splits && + most_splits <= SL_MAX_SPLITS / 2) { + col_price = prices->get(prices->context); + } + if (col_price) { + cost lb = 0; + v_init.clear(); + for (split_int bi = 0; bi < b.n_splits; ++bi) { + if (!b_match[bi]) { + lb += col_price[bi]; + v_init.push_back(col_price[bi]); + } + } + for (split_int ai = 0; ai < a.n_splits; ++ai) { + if (a_match[ai]) continue; + const cost *row = score.row(ai); + cost row_min = row[0] - col_price[0]; + for (split_int bi = 1; bi < b.n_splits; ++bi) { + const cost reduced = row[bi] - col_price[bi]; + if (reduced < row_min) row_min = reduced; + } + lb += row_min; + } + const cost lap_ub = max_score * lap_dim - lb; + const double bound = exact_matches ? + static_cast(lap_ub) * over_max_score + + exact_match_score / n_tips : + static_cast(lap_ub) / max_score; + if (bound < min_score) { + return List::create(Named("score") = bound, + _["matching"] = IntegerVector(0)); + } + } + } + std::vector rowsol; std::vector colsol; resize_uninitialized(rowsol, lap_dim); resize_uninitialized(colsol, lap_dim); + [[maybe_unused]] static thread_local LapScratch warm_scratch; + const auto Solve = [&](cost_matrix &m) -> cost { + if constexpr (Bounded) { + if (col_price) { + return lap_warm(lap_dim, m, rowsol, colsol, warm_scratch, + v_init.data()); + } + } + return lap(lap_dim, m, rowsol, colsol, !Bounded); + }; cost_matrix local_small(Bounded ? 0 : lap_dim); cost_matrix &small_score = Bounded ? BoundedPool(1, lap_dim) : local_small; @@ -604,7 +676,7 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, small_score.padAfterRow(lap_dim - b_extra_splits, max_score); const double lap_score = static_cast((max_score * lap_dim) - - lap(lap_dim, small_score, rowsol, colsol, !Bounded)) * over_max_score; + Solve(small_score)) * over_max_score; const double final_score = lap_score + (exact_match_score / n_tips); std::unique_ptr lap_decode = std::make_unique(lap_dim); @@ -642,7 +714,7 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, } const double final_score = static_cast( - (max_score * lap_dim) - lap(lap_dim, score, rowsol, colsol, !Bounded) + (max_score * lap_dim) - Solve(score) ) / max_score; std::vector final_matching; @@ -661,9 +733,11 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, template List mutual_clustering(const RawMatrix &x, const RawMatrix &y, - const int32 n_tips, const double min_score = 0) { + const int32 n_tips, const double min_score = 0, + const ColPrices *prices = nullptr) { return TreeDist::dispatch_popcount([&](auto pop) { - return mutual_clustering_impl(x, y, n_tips, min_score); + return mutual_clustering_impl(x, y, n_tips, + min_score, prices); }); } @@ -822,9 +896,71 @@ List cpp_mutual_clustering(const RawMatrix &x, const RawMatrix &y, // min_score. Used to discard candidates early in SPRDist(method = "Rogue"). double mutual_clustering_bounded_score(const RawMatrix &x, const RawMatrix &y, const int32 n_tip, - const double min_score) { + const double min_score, + const ColPrices *prices) { return Rcpp::as( - mutual_clustering(x, y, n_tip, min_score)["score"]); + mutual_clustering(x, y, n_tip, min_score, prices)["score"]); +} + +// Optimal column duals of the mutual clustering LAP between the splits `x` +// and `y` (identical splits included), one per split in `y`, in cost units, +// shifted so that the least is zero and capped at BIG: prices for +// mutual_clustering() on related problems. Empty if there is nothing +// to match. +std::vector mutual_clustering_col_prices(const RawMatrix &x, + const RawMatrix &y, + const int32 n_tips) { + const SplitList a(x); + const SplitList b(y); + const split_int dim = std::max(a.n_splits, b.n_splits); + if (dim == 0 || n_tips < 4 || b.n_splits == 0) { + return {}; // #nocov + } + constexpr cost max_score = BIG; + const double max_over_tips = static_cast(max_score) / n_tips; + const double lg2_n = lg2_lookup(n_tips); + cost_matrix score(dim); + for (split_int ai = 0; ai < a.n_splits; ++ai) { + const split_int na = a.in_split[ai]; + const split_int nA = n_tips - na; + const double offset_a = lg2_n - lg2_lookup(na); + const double offset_A = lg2_n - lg2_lookup(nA); + for (split_int bi = 0; bi < b.n_splits; ++bi) { + split_int a_and_b = 0; + for (split_int bin = 0; bin < a.n_bins; ++bin) { + a_and_b += TreeDist::count_bits(a.state[ai][bin] & b.state[bi][bin]); + } + const split_int nb = b.in_split[bi]; + const split_int nB = n_tips - nb; + const split_int a_and_B = na - a_and_b; + const split_int A_and_b = nb - a_and_b; + const split_int A_and_B = nA - A_and_b; + if (a_and_b == A_and_b && a_and_b == a_and_B && a_and_b == A_and_B) { + score(ai, bi) = max_score; + } else { + const double lg2_nb = lg2_lookup(nb); + const double lg2_nB = lg2_lookup(nB); + const double ic_sum = + a_and_b * (lg2_lookup(a_and_b) + offset_a - lg2_nb) + + a_and_B * (lg2_lookup(a_and_B) + offset_a - lg2_nB) + + A_and_b * (lg2_lookup(A_and_b) + offset_A - lg2_nb) + + A_and_B * (lg2_lookup(A_and_B) + offset_A - lg2_nB); + score(ai, bi) = max_score - static_cast(ic_sum * max_over_tips); + } + } + score.padRowAfterCol(ai, b.n_splits, max_score); + } + score.padAfterRow(a.n_splits, max_score); + std::vector rowsol(dim); + std::vector colsol(dim); + LapScratch scratch; + lap(dim, score, rowsol, colsol, false, scratch); + std::vector ret(scratch.v.begin(), scratch.v.begin() + b.n_splits); + const cost least = *std::min_element(ret.begin(), ret.end()); + for (cost &v : ret) { + v = std::min(v - least, max_score); + } + return ret; } // [[Rcpp::export]] @@ -834,7 +970,7 @@ double cpp_mutual_clustering_bounded(const RawMatrix &x, const RawMatrix &y, ASSERT(x.cols() == y.cols() && "Input splits must address same number of tips."); const int32 n_tip = static_cast(nTip[0]); TreeDist::check_ntip(n_tip); - return mutual_clustering_bounded_score(x, y, n_tip, minScore); + return mutual_clustering_bounded_score(x, y, n_tip, minScore, nullptr); } // [[Rcpp::export]]