fix(stats): aggregate DuckDB counts exactly

.github/workflows/ci.yml

@@ -400,11 +400,11 @@ 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/validated-statistics-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/validated-statistics-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}'; }
-            echo "::error title=Validated-statistics tests failed::$(tail -n 120 /tmp/validated-statistics-tests.log | enc)"
+            echo "::error title=Validated-statistics tests failed::$(tail -c 3500 /tmp/validated-statistics-tests.log | enc)"
           fi
           exit "$status"
 

src/analysis/clade_cumulus.jl

@@ -28,7 +28,9 @@ const DEFAULT_TAXONOMIC_RANKS = [
     "Class", "Order", "Family", "Genus", "Species",
 ]
 
-const PLACEMENT_WEIGHT_SUM_TOL = 1e-5
+# Placement weights are often rounded or truncated in serialized jplace output.
+# Accept small omissions without treating materially incomplete evidence as valid.
+const PLACEMENT_WEIGHT_SUM_TOL = 0.05
 
 ## Phylogenetic tree structures & Newick / jplace parsers
 mutable struct PhyloNode
@@ -202,8 +204,11 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]):
         end
     end
 
-    function read_label_and_quoted()::Tuple{String,Bool}
-        skip_ws_and_leading_comments!()
+    function read_label_and_quoted(;
+        capture_nhx::Union{Ref{Union{Float64,Nothing}},Nothing} = nothing,
+        capture_edge::Union{Ref{Union{Int,Nothing}},Nothing} = nothing,
+    )::Tuple{String,Bool}
+        skip_ws_and_leading_comments!(capture_nhx = capture_nhx, capture_edge = capture_edge)
         pos[] > n && return ("", false)
         if bytes[pos[]] == UInt8('\'')
             pos[] += 1
@@ -260,7 +265,8 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]):
 
         nhx_sup = Ref{Union{Float64,Nothing}}(nothing)
         edge_num = Ref{Union{Int,Nothing}}(nothing)
-        raw_lbl, was_quoted = read_label_and_quoted()
+        raw_lbl, was_quoted = read_label_and_quoted(
+            capture_nhx = nhx_sup, capture_edge = edge_num)
         skip_ws_and_leading_comments!(capture_nhx = nhx_sup, capture_edge = edge_num)
 
         branch_len = 0.0
@@ -991,7 +997,7 @@ function _graft_jplace_queries(
 
     for q in jd.queries
         q.status === :valid && !isempty(q.candidates) || continue
-        total_w = sum((c.weight for c in q.candidates); init = 0.0)
+        total_w = sum((c.like_weight for c in q.candidates); init = 0.0)
         total_w > 0.0 || continue
         multi = length(q.candidates) > 1
         for (qname, mult) in q.names
@@ -1000,7 +1006,7 @@ function _graft_jplace_queries(
             tx = get(asv_deepest_taxon, qname, nothing)
             tconf = get(asv_tax_conf, qname, 1.0)
             for (ci, cand) in enumerate(q.candidates)
-                share = cand.weight / total_w
+                share = cand.like_weight / total_w
                 tip_id = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname
                 if multi
                     asv_reads[tip_id] = get(asv_reads, tip_id, 0.0) + base_cnt * share
@@ -1009,9 +1015,10 @@ function _graft_jplace_queries(
                 elseif !haskey(asv_reads, tip_id)
                     asv_reads[tip_id] = base_cnt
                 end
-                query_lwr[tip_id] = cand.weight
+                query_lwr[tip_id] = cand.like_weight
                 push!(get!(by_edge, cand.edge_num, _GraftedQuerySpec[]),
-                    _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_length), max(0.0, cand.pendant_length), cand.weight))
+                    _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_len),
+                        max(0.0, cand.pendant_len), cand.like_weight))
             end
         end
     end
@@ -1021,8 +1028,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.
@@ -1040,7 +1047,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)
@@ -1815,6 +1822,9 @@ function _reaggregate_after_move(
                 merged_taxa[k] = get(merged_taxa, k, 0) + v
             end
         end
+        # Match clade_faith_pd: an active non-root clade includes the branch
+        # connecting it to its parent as well as its active descendant branches.
+        any_active && !is_root && (sub_pd += max(0.0, u.branch_length))
 
         dom_tx = nothing
         dom_n = 0

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

@@ -93,8 +93,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,
@@ -200,7 +200,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"))))
@@ -268,13 +268,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
@@ -391,7 +391,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"))))

src/server/routes/validated_statistics.jl

@@ -38,18 +38,71 @@ function _statistics_guard(f::Function)
     end
 end
 
-function _statistics_body(req)
+function _statistics_body(req; reject_noninteger_tokens::Bool=false)
     length(req.body) <= _STATISTICS_REQUEST_LIMIT || throw(ResourceLimitError(
         "request bytes", _STATISTICS_REQUEST_LIMIT, length(req.body)))
+    source = String(req.body)
     parsed = try
-        JSON3.read(String(req.body))
+        JSON3.read(source)
     catch
         throw(ArgumentError("request body is not valid JSON"))
     end
     parsed isa JSON3.Object || throw(ArgumentError("request body must be a JSON object"))
+    reject_noninteger_tokens && _reject_noninteger_json_tokens(source)
     parsed
 end
 
+"""Reject decimal/exponent tokens whose representation JSON3 would erase.
+JSON3 1.14 normalises integral `1.0` and `1e0` tokens to Int64. Every numeric
+field in a direct validated-statistics request is contractually integral, so a
+lexical pass over the already-validated JSON preserves that boundary."""
+function _reject_noninteger_json_tokens(source::AbstractString)
+    bytes = codeunits(source)
+    position = 1
+    while position <= length(bytes)
+        byte = bytes[position]
+        if byte == UInt8('"')
+            position += 1
+            closed = false
+            while position <= length(bytes)
+                byte = bytes[position]
+                if byte == UInt8('\\')
+                    position += 2 # JSON3 has already validated the escape.
+                elseif byte == UInt8('"')
+                    closed = true
+                    position += 1
+                    break
+                else
+                    position += 1
+                end
+            end
+            closed || throw(ArgumentError("unterminated JSON string"))
+        elseif 0x30 <= byte <= 0x39 || byte == UInt8('-') || byte == UInt8('+')
+            token_start = position
+            noninteger = false
+            while position <= length(bytes)
+                byte = bytes[position]
+                if 0x30 <= byte <= 0x39 || byte == UInt8('+') || byte == UInt8('-') ||
+                   byte == UInt8('.') || byte == UInt8('e') || byte == UInt8('E')
+                    noninteger |= byte == UInt8('.') || byte == UInt8('e') || byte == UInt8('E')
+                    position += 1
+                else
+                    break
+                end
+            end
+            if noninteger
+                token = String(bytes[token_start:position - 1])
+                throw(UnsupportedRepresentationError(
+                    "an integer JSON token or whole-number string", token,
+                    "validated direct statistics; decimal and exponent JSON numbers are not exact input"))
+            end
+        else
+            position += 1
+        end
+    end
+    nothing
+end
+
 function _policy_integer(raw, key::Symbol, default)
     value = get(raw, key, default)
     (value isa Integer && !(value isa Bool)) || throw(ArgumentError(
@@ -95,8 +148,8 @@ function _json_exact_count(value, feature::String, sample::String)
                 "feature '$feature', sample '$sample'; use a canonical whole-number string"))
         end
     end
-    # JSON3 represents a decimal/exponent token as a float even when its value
-    # is mathematically integral. Reject the representation, not just 1.5.
+    # Non-integral decimal/exponent tokens remain floats after JSON3 parsing.
+    # Integral ones are rejected by the raw-token pass before this function.
     (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError(
         "an integer JSON token or whole-number string", repr(value),
         "feature '$feature', sample '$sample'; floating-point JSON is not exact input"))
@@ -180,8 +233,6 @@ function _validated_integer_sample_columns(con, table::String, columns::Vector{S
     samples
 end
 
-"""Read and aggregate an integer-safe run table. Only GROUP BY is done in SQL;
-all count addition is Julia BigInt. No SQL SUM and no Float64 matrix is used."""
 function _reject_float_count_filters(sample_columns::Vector{String}, params::Dict{String,String})
     for sample in sample_columns, prefix in ("col_min.", "col_max.")
         key = prefix * sample
@@ -192,6 +243,9 @@ function _reject_float_count_filters(sample_columns::Vector{String}, params::Dic
     nothing
 end
 
+"""Read and aggregate an integer-safe run table. SQL performs only projection
+and filtering; feature grouping and all count addition use Julia BigInt. No SQL
+SUM, arbitrary grouped value, or Float64 matrix is used."""
 function _validated_duckdb_counts(con, table::String, columns::Vector{String};
                                   prefix=nothing, taxonomic_level=nothing,
                                   where_clause::String="", where_params=Any[])
@@ -200,6 +254,7 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String};
         samples = filter(c -> startswith(c, prefix), samples)
     end
     isempty(samples) && throw(ArgumentError("no sample columns remain after the prefix filter"))
+    sample_symbols = Symbol.(samples)
 
     feature_expression = if isnothing(taxonomic_level)
         identity = Analysis.sequence_column_name(columns)
@@ -219,29 +274,67 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String};
     end
     count_columns = join([_sql_identifier(c) for c in samples], ", ")
     sql = "SELECT $feature_expression AS feature, $count_columns FROM $(_sql_identifier(table)) " *
-          "$where_clause GROUP BY $feature_expression ORDER BY feature LIMIT $(MAX_CELLS + 1)"
-    frame = DataFrame(DBInterface.execute(con, sql, where_params))
-    nrow(frame) <= MAX_CELLS || throw(ResourceLimitError("features", MAX_CELLS, nrow(frame)))
-    nrow(frame) > 0 || throw(ArgumentError("no rows remain after filtering"))
-    nrow(frame) * length(samples) <= MAX_CELLS || throw(ResourceLimitError(
-        "table cells", MAX_CELLS, nrow(frame) * length(samples)))
+          where_clause
 
-    features = Vector{String}(undef, nrow(frame))
-    counts = Matrix{BigInt}(undef, nrow(frame), length(samples))
-    for (i, row) in enumerate(eachrow(frame))
-        ismissing(row.feature) && throw(UnsupportedRepresentationError(
-            "a non-missing feature label", "missing", "validated run statistics row $i"))
-        features[i] = string(row.feature)
-        for (j, sample) in enumerate(samples)
-            value = row[Symbol(sample)]
-            ismissing(value) && throw(UnsupportedRepresentationError(
-                "a non-missing integer count", "NULL", "feature '$(features[i])', sample '$sample'"))
-            (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError(
-                "an integer count", repr(value), "feature '$(features[i])', sample '$sample'"))
-            value >= 0 || throw(ArgumentError(
-                "negative count $value at feature '$(features[i])', sample '$sample'"))
-            counts[i, j] = BigInt(value) # every addition later is exact BigInt
+    # A SQL GROUP BY would either overflow fixed-width SUM or require an
+    # arbitrary per-group value. Stream the raw rows instead, and bound the
+    # aggregated matrix as features arrive rather than materialising the raw
+    # table in Julia. DuckDB.Tables is used through the direct DuckDB dependency.
+    features = String[]
+    rows_by_feature = Vector{Vector{BigInt}}()
+    feature_rows = Dict{String,Int}()
+    raw_rows = 0
+    statement = DBInterface.prepare(con, sql, DuckDB.StreamResult)
+    try
+        result = DBInterface.execute(statement, where_params)
+        try
+            for chunk in DuckDB.Tables.partitions(result)
+                chunk_columns = DuckDB.Tables.columns(chunk)
+                for row in DuckDB.Tables.rows(chunk_columns)
+                    raw_rows += 1
+                    ismissing(row.feature) && throw(UnsupportedRepresentationError(
+                        "a non-missing feature label", "missing",
+                        "validated run statistics row $raw_rows"))
+                    feature = string(row.feature)
+                    i = get(feature_rows, feature, 0)
+                    if iszero(i)
+                        n_features = length(features) + 1
+                        n_features <= MAX_CELLS || throw(ResourceLimitError(
+                            "features", MAX_CELLS, n_features))
+                        n_features * length(samples) <= MAX_CELLS || throw(ResourceLimitError(
+                            "table cells", MAX_CELLS, n_features * length(samples)))
+                        push!(features, feature)
+                        push!(rows_by_feature, zeros(BigInt, length(samples)))
+                        i = length(features)
+                        feature_rows[feature] = i
+                    end
+                    for (j, sample) in enumerate(samples)
+                        value = row[sample_symbols[j]]
+                        ismissing(value) && throw(UnsupportedRepresentationError(
+                            "a non-missing integer count", "NULL",
+                            "feature '$feature', sample '$sample'"))
+                        (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError(
+                            "an integer count", repr(value), "feature '$feature', sample '$sample'"))
+                        value >= 0 || throw(ArgumentError(
+                            "negative count $value at feature '$feature', sample '$sample'"))
+                        rows_by_feature[i][j] += BigInt(value)
+                    end
+                end
+            end
+        finally
+            DBInterface.close!(result)
         end
+    finally
+        DBInterface.close!(statement)
+    end
+
+    isempty(features) && throw(ArgumentError("no rows remain after filtering"))
+    order = sortperm(features)
+    features = features[order]
+    rows_by_feature = rows_by_feature[order]
+    counts = Matrix{BigInt}(undef, length(features), length(samples))
+    for i in eachindex(features), j in eachindex(samples)
+        counts[i, j] = rows_by_feature[i][j]
     end
     ExactSummaries.validate_coordinates(counts, samples, features)
     counts, samples, features
@@ -280,7 +373,7 @@ end
 
 @post "/api/v1/statistics/summaries" function(req)
     _statistics_guard() do
-        body = _statistics_body(req)
+        body = _statistics_body(req; reject_noninteger_tokens=true)
         policy = _statistics_policy(body)
         counts, labels, features, groups = _direct_count_input(body)
         result = _run_statistics(counts, labels, features, groups, policy)