Land #217 (exact DuckDB counts) on main, and fix the Wald tail

.github/workflows/ci.yml

@@ -400,7 +400,7 @@ jobs:
         shell: bash
         run: |
           set +e
-          julia --project=. -t2 --compiled-modules=no test/runtests.jl numeric_policy exact_summaries high_precision derived_analyses validated_statistics 2>&1 | tee /tmp/julia-validated-tests.log
+          julia --project=. -t2 --compiled-modules=no test/runtests.jl validated_statistics numeric_policy exact_summaries high_precision derived_analyses 2>&1 | tee /tmp/julia-validated-tests.log
           status="${PIPESTATUS[0]}"
           if [ "$status" -ne 0 ]; then
             enc() { head -c 3500 | sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; }

src/analysis/clade_cumulus.jl

@@ -1029,8 +1029,8 @@ function _graft_jplace_queries(
 
     function walk(u::PhyloNode)::PhyloNode
         rebuilt_kids = PhyloNode[walk(ch) for ch in u.children]
-        cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_num, rebuilt_kids)
-        specs = u.edge_num === nothing ? nothing : get(by_edge, u.edge_num, nothing)
+        cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_number, rebuilt_kids)
+        specs = u.edge_number === nothing ? nothing : get(by_edge, u.edge_number, nothing)
         (specs === nothing || isempty(specs)) && return cur
 
         # Order grafts from most distal (closest to `cur`, smallest `distal`) to most proximal.
@@ -1048,7 +1048,7 @@ function _graft_jplace_queries(
             sub = PhyloNode(string("jplace_clade#", seq_auto[]), "", 0.0, 1.0, nothing, PhyloNode[sub, qtip])
         end
         rem_len = max(0.0, total_len - prev_distal)
-        PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_num, sub.children)
+        PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_number, sub.children)
     end
 
     walk(base_root)
@@ -1791,11 +1791,15 @@ function _reaggregate_after_move(
                 mono = get(total_tips_per_taxon, u.dominant_taxon, 0) == 1
             end
             active_here = u.cumulative_count > 0.0
+            # `subtree_pd` is node-local, as in `clade_faith_pd`: a node's own
+            # (non-root) branch counts once, here, and its parent adds only the
+            # value it returns.
+            own_pd = (active_here && !is_root) ? max(0.0, u.branch_length) : 0.0
             st = Epistemic.epistemic_status(case_obj, u.id)
             leaf_node = CladeNode(
                 u.id, u.name, u.rank, u.direct_count, u.cumulative_count, freq,
                 u.residual_count, Epistemic.cloud_size_by_residual(u.residual_count),
-                u.branch_length, u.support, 0.0, u.dominant_taxon, u.taxon_purity, mono,
+                u.branch_length, u.support, own_pd, u.dominant_taxon, u.taxon_purity, mono,
                 st, Epistemic.epistemic_colour(st), CladeNode[],
             )
             return (leaf_node, tmap, 1, active_here)
@@ -1817,12 +1821,18 @@ function _reaggregate_after_move(
             n_sub_tips += nt
             if ch_active
                 any_active = true
-                sub_pd += max(0.0, rk.branch_length) + rk.subtree_pd
+                # The child's own branch is already inside `rk.subtree_pd`;
+                # adding `rk.branch_length` again would count every active
+                # internal edge twice.
+                sub_pd += rk.subtree_pd
             end
             for (k, v) in ct
                 merged_taxa[k] = get(merged_taxa, k, 0) + v
             end
         end
+        # Match clade_faith_pd: an active non-root clade adds the branch
+        # connecting it to its parent, once, on top of its descendants' values.
+        any_active && !is_root && (sub_pd += max(0.0, u.branch_length))
 
         dom_tx = nothing
         dom_n = 0

src/analysis/compositional.jl

@@ -42,7 +42,8 @@
 # SpecialFunctions, so the models add no dependency: `_loggamma` and `_digamma`
 # are the Stirling series after the recurrence is applied up to x >= 10,
 # accurate to about 1e-14 (checked against Julia's `lgamma` and known constants
-# in the unit tests).
+# in the unit tests), and `_erfc` is the regularised incomplete gamma the Wald
+# p-value needs.
 module CompositionalModels
 
 export CompositionalFit, CompositionalRefusal, RuntimeRefusal, MAX_TAXA,
@@ -193,6 +194,57 @@ function compositional_digamma(x::Real)
         inv2 * (1 / 252 - inv2 * (1 / 240 - inv2 / 132)))) - shift
 end
 
+# The Wald test reports the two-sided normal tail, `erfc(|z| / sqrt(2))`. The
+# complementary error function is the regularised upper incomplete gamma
+# Q(1/2, t^2), so the power series of P(1/2, x) evaluates it below the switch
+# and the continued fraction of Q(1/2, x) above it (modified Lentz). Each branch
+# ends within a few dozen terms at the switch and fewer away from it, and the
+# result is accurate to about 1e-14 relative over the range the Wald statistic
+# reaches (verified against independent reference values in the unit tests).
+const _LOG_SQRT_PI = 0.5723649429247000870717136756765293558
+
+function _erfc(t::Float64)::Float64
+    t < 0.0 && return 2.0 - _erfc(-t)
+    iszero(t) && return 1.0
+    x = t * t
+    # exp(-x) * x^(1/2) / Gamma(1/2); x can underflow to zero for tiny t, and
+    # then exp(-Inf) is zero and the answer is 1 to double precision anyway.
+    log_prefactor = -x + 0.5 * log(x) - _LOG_SQRT_PI
+    if x < 1.5
+        # P(1/2, x) as a power series, then erfc = 1 - P. Below the switch the
+        # subtraction is safe: P <= 0.92 and erfc >= 0.08.
+        ap = 0.5
+        term = 1.0 / ap
+        acc = term
+        for _ in 1:200
+            ap += 1.0
+            term *= x / ap
+            acc += term
+            abs(term) <= abs(acc) * 1e-16 && break
+        end
+        return 1.0 - acc * exp(log_prefactor)
+    end
+    # Q(1/2, x) by its continued fraction.
+    tiny = 1e-300
+    b = x + 0.5
+    c = 1.0 / tiny
+    d = 1.0 / b
+    h = d
+    for i in 1:200
+        an = -float(i) * (float(i) - 0.5)
+        b += 2.0
+        d = an * d + b
+        abs(d) < tiny && (d = tiny)
+        c = b + an / c
+        abs(c) < tiny && (c = tiny)
+        d = 1.0 / d
+        delta = d * c
+        h *= delta
+        abs(delta - 1.0) <= 1e-16 && break
+    end
+    exp(log_prefactor) * h
+end
+
 ## Geometry
 
 """
@@ -860,7 +912,7 @@ function _wald(coefficients::Matrix{Float64}, variances::Vector{Float64})
         se[l, k] = s
         if s > 0.0 && isfinite(s)
             z[l, k] = coefficients[l, k] / s
-            p[l, k] = erfc(abs(z[l, k]) / sqrt(2.0))
+            p[l, k] = _erfc(abs(z[l, k]) / sqrt(2.0))
         else
             z[l, k] = NaN
             p[l, k] = NaN

src/pipeline/phylogeny.jl

@@ -33,6 +33,7 @@ export REFERENCE_STEPS, PLACEMENT_STEPS, Step, placement_settings, run_workflow,
     using SHA, JSON3, Dates, Logging, OrderedCollections
     using ..PipelineLog
     using ..Validation
+    import ..Tools
     using ..Tools: tool_bin, _sq, _run_logged, _safe_optional_args, _num
     using ..RemoteExec: run_remote, checked_target
     import ..Provenance

src/server/routes/clade_cumulus.jl

@@ -97,8 +97,8 @@ function _resolve_tree_inputs_from_request(
     source_digests::OrderedDict{String,String},
 )
     tree_source = lowercase(strip(string(get(body, :treeSource, "auto"))))
-    tree_file = strip(string(get(body, :treeFile, "")))
-    inline_newick = strip(string(get(body, :newick, "")))
+    tree_file = String(strip(string(get(body, :treeFile, ""))))
+    inline_newick = String(strip(string(get(body, :newick, ""))))
 
     if tree_source == "study_tree" || (!isempty(tree_file) && tree_source == "auto")
         _valid_tree_file(tree_file) || return (nothing, nothing, :auto,
@@ -208,7 +208,7 @@ end
         return json_error(400, "invalid_json", "Request body must be valid JSON")
     end
 
-    table = strip(string(get(body, :table, "merged")))
+    table = String(strip(string(get(body, :table, "merged"))))
     _valid_name(table) || return json_error(400, "invalid_table", "Invalid table name '$table'")
 
     mode_str = lowercase(strip(string(get(body, :mode, "combined"))))
@@ -276,13 +276,13 @@ end
             jd_div = MetaManifold.CladeCumulus.parse_jplace(jplace_str)
             for q in jd_div.queries
                 q.status === :valid && !isempty(q.candidates) || continue
-                tw = sum((c.weight for c in q.candidates); init = 0.0)
+                tw = sum((c.like_weight for c in q.candidates); init = 0.0)
                 tw > 0.0 || continue
                 multi = length(q.candidates) > 1
                 for (qname, _) in q.names
                     for (ci, cand) in enumerate(q.candidates)
                         tid = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname
-                        tip_weights[tid] = cand.weight / tw
+                        tip_weights[tid] = cand.like_weight / tw
                     end
                 end
             end
@@ -399,7 +399,7 @@ end
     (isempty(source_id) || isempty(target_id)) &&
         return json_error(400, "missing_clade_ids", "Both sourceId and targetId are required")
 
-    table = strip(string(get(body, :table, "merged")))
+    table = String(strip(string(get(body, :table, "merged"))))
     _valid_name(table) || return json_error(400, "invalid_table", "Invalid table name '$table'")
 
     mode_str = lowercase(strip(string(get(body, :mode, "combined"))))