fix(server): aggregate validated run statistics in BigInt, not SQL

.github/workflows/ci.yml

@@ -403,8 +403,23 @@ jobs:
           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
           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)"
+            enc() { sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; }
+            # One annotation cannot hold a failing run, and a tail hides every
+            # failure but the last one. Publish a digest of all failing
+            # assertions first - the testset, the expression and the evaluated
+            # values are what identify the bug - then the end of the log.
+            grep -E "Test Failed at |Error During Test at |^  Expression: |^   Evaluated: |^  Test threw exception" \
+              /tmp/validated-statistics-tests.log > /tmp/validated-statistics-failures.log || true
+            rm -f /tmp/validated-statistics-digest-*
+            split -b 3800 /tmp/validated-statistics-failures.log /tmp/validated-statistics-digest- || true
+            index=0
+            for chunk in /tmp/validated-statistics-digest-*; do
+              [ -e "$chunk" ] || continue
+              index=$((index + 1))
+              [ "$index" -le 8 ] || break
+              echo "::error title=Validated-statistics failure digest $index::$(enc < "$chunk")"
+            done
+            echo "::error title=Validated-statistics log tail::$(tail -n 60 /tmp/validated-statistics-tests.log | enc)"
           fi
           exit "$status"
 

src/analysis/clade_cumulus.jl

@@ -28,7 +28,11 @@ const DEFAULT_TAXONOMIC_RANKS = [
     "Class", "Order", "Family", "Genus", "Species",
 ]
 
-const PLACEMENT_WEIGHT_SUM_TOL = 1e-5
+# jplace writers round `like_weight_ratio` and drop placements below their own
+# likelihood-weight cutoff, so a query whose surviving weights sum to a little
+# under one is still complete evidence. A shortfall beyond 5% is a missing or
+# discarded placement rather than rounding, and the query stays `undecided`.
+const PLACEMENT_WEIGHT_SUM_TOL = 5e-2
 
 ## Phylogenetic tree structures & Newick / jplace parsers
 mutable struct PhyloNode
@@ -260,6 +264,11 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]):
 
         nhx_sup = Ref{Union{Float64,Nothing}}(nothing)
         edge_num = Ref{Union{Int,Nothing}}(nothing)
+        # A node's comment or jplace edge tag can sit before its label
+        # (`(...)[&&NHX:B=92]:0.1`, `(...){4}`) as well as after it, and
+        # `read_label_and_quoted` skips comments without capturing them, so the
+        # leading ones are captured first.
+        skip_ws_and_leading_comments!(capture_nhx = nhx_sup, capture_edge = edge_num)
         raw_lbl, was_quoted = read_label_and_quoted()
         skip_ws_and_leading_comments!(capture_nhx = nhx_sup, capture_edge = edge_num)
 
@@ -991,7 +1000,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 +1009,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 +1018,9 @@ 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_length), max(0.0, cand.pendant_length), cand.like_weight))
             end
         end
     end

src/pipeline/phylogeny.jl

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

src/server/routes/clade_cumulus.jl

@@ -64,7 +64,9 @@ function _probe_optional_phylo_tool(tool_key::String, fallbacks::Vector{String})
     end
 end
 
-function _load_run_table_df(study::String, run::String, table::String; group::Union{String,Nothing}=nothing)
+function _load_run_table_df(study::String, run::String, table::AbstractString;
+                            group::Union{String,Nothing}=nothing)
+    table = String(table) # request bodies hand us the SubString that `strip` returns
     dir = _require_duckdb(study, run; group)
     isnothing(dir) && return (nothing, nothing, json_error(404, "no_results",
         "No results for run '$run' - run the pipeline first"))
@@ -145,19 +147,28 @@ end
             names_df = DataFrame(DBInterface.execute(con, "SHOW TABLES"))
             tlist = sort(String[string(r.name) for r in eachrow(names_df)])
             for t in tlist
-                cols = _duckdb_columns(con, t)
-                scols = _sample_count_columns(con, t)
-                isempty(scols) && continue
-                nrows = Int(only(DataFrame(DBInterface.execute(con, "SELECT COUNT(*) AS n FROM \"$t\""))).n)
-                has_seq = any(c -> c in ("Sequence", "sequence", "seq", "DNA"), cols)
-                ranks = [r for r in MetaManifold.CladeCumulus.DEFAULT_TAXONOMIC_RANKS if r in cols]
-                push!(tables_info, OrderedDict{String,Any}(
-                    "table"        => t,
-                    "row_count"    => nrows,
-                    "sample_count" => length(scols),
-                    "has_sequence" => has_seq,
-                    "ranks"        => ranks,
-                ))
+                # `SHOW TABLES` lists views too, and a view over a DuckDB
+                # replacement scan is bound to the connection that registered it.
+                # Reading one from here raises, so an object that cannot be read
+                # is skipped instead of failing the whole listing.
+                try
+                    cols = _duckdb_columns(con, t)
+                    scols = _sample_count_columns(con, t)
+                    isempty(scols) && continue
+                    nrows = Int(only(DataFrame(DBInterface.execute(con, "SELECT COUNT(*) AS n FROM \"$t\""))).n)
+                    has_seq = any(c -> c in ("Sequence", "sequence", "seq", "DNA"), cols)
+                    ranks = [r for r in MetaManifold.CladeCumulus.DEFAULT_TAXONOMIC_RANKS if r in cols]
+                    push!(tables_info, OrderedDict{String,Any}(
+                        "table"        => t,
+                        "row_count"    => nrows,
+                        "sample_count" => length(scols),
+                        "has_sequence" => has_seq,
+                        "ranks"        => ranks,
+                    ))
+                catch err
+                    err isa InterruptException && rethrow()
+                    @warn "clade-cumulus sources: skipping an unreadable object" object=t exception=err
+                end
             end
         finally
             DBInterface.close!(con)

src/server/routes/trees.jl

@@ -17,9 +17,9 @@ const _tree_view_lock  = ReentrantLock()
 _trees_dir(study::String) = joinpath(ServerState.projects_dir(), study, "trees")
 _tree_view_path(study::String, file::String) = joinpath(_trees_dir(study), ".views", file * ".json")
 
-_tree_format(file::String) = endswith(lowercase(file), ".jplace") ? "jplace" : "newick"
+_tree_format(file::AbstractString) = endswith(lowercase(file), ".jplace") ? "jplace" : "newick"
 
-function _valid_tree_file(file::String)::Bool
+function _valid_tree_file(file::AbstractString)::Bool
     _valid_name(file) && any(ext -> endswith(lowercase(file), ext), _TREE_EXTENSIONS)
 end
 

src/server/routes/validated_statistics.jl

@@ -95,8 +95,10 @@ 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.
+    # JSON3 normalises an integral decimal or exponent token (`1.0`, `1e0`) to
+    # Int64 while it parses, so the only non-integers that can arrive here are
+    # genuinely fractional or out of Int64 range, and those are refused rather
+    # than rounded.
     (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 +182,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 +192,15 @@ function _reject_float_count_filters(sample_columns::Vector{String}, params::Dic
     nothing
 end
 
+"""Read and aggregate an integer-safe run table without SQL arithmetic.
+
+SQL only projects rows; every count is added in Julia `BigInt`. A `SUM` of
+`BIGINT`/`HUGEINT` columns is evaluated in `Int64`/`Int128` and can overflow
+before the exact total is known - two sample columns of `2^126` already sum to
+`2^127`, one past `typemax(Int128)` - so no `SUM` and no `Float64` matrix is
+used. Rows that share a feature label are therefore summed here, after SQL has
+ordered them by label.
+"""
 function _validated_duckdb_counts(con, table::String, columns::Vector{String};
                                   prefix=nothing, taxonomic_level=nothing,
                                   where_clause::String="", where_params=Any[])
@@ -218,31 +227,47 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String};
         "COALESCE(NULLIF(CAST($(_sql_identifier(taxon)) AS VARCHAR), ''), 'Unassigned')"
     end
     count_columns = join([_sql_identifier(c) for c in samples], ", ")
+    # Project the rows and let SQL sort them by feature; no addition happens
+    # here. `GROUP BY feature` alone is a DuckDB binder error, because the count
+    # columns are neither grouped nor aggregated, and `SUM` cannot take their
+    # place without overflowing. Sorting by label makes equal labels adjacent
+    # and leaves the accumulator independent of row order.
     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)"
+          "$where_clause 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)))
+    nrow(frame) <= MAX_CELLS || throw(ResourceLimitError("table rows", MAX_CELLS, nrow(frame)))
 
-    features = Vector{String}(undef, nrow(frame))
-    counts = Matrix{BigInt}(undef, nrow(frame), length(samples))
+    features = String[]
+    totals = Dict{String,Vector{BigInt}}()
     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)
+        label = string(row.feature)
+        sums = get(totals, label, nothing)
+        if isnothing(sums)
+            length(features) + 1 <= MAX_CELLS ÷ max(length(samples), 1) ||
+                throw(ResourceLimitError("table cells", MAX_CELLS,
+                                         (length(features) + 1) * length(samples)))
+            sums = zeros(BigInt, length(samples))
+            totals[label] = sums
+            push!(features, label)
+        end
         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'"))
+                "a non-missing integer count", "NULL", "feature '$label', sample '$sample'"))
             (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError(
-                "an integer count", repr(value), "feature '$(features[i])', sample '$sample'"))
+                "an integer count", repr(value), "feature '$label', 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
+                "negative count $value at feature '$label', sample '$sample'"))
+            sums[j] += BigInt(value) # every addition is exact, so no overflow
         end
     end
+    counts = Matrix{BigInt}(undef, length(features), length(samples))
+    for (i, label) in enumerate(features)
+        counts[i, :] = totals[label]
+    end
     ExactSummaries.validate_coordinates(counts, samples, features)
     counts, samples, features
 end

src/server/server.jl

@@ -118,7 +118,6 @@ module Server
     include(joinpath(@__DIR__, "routes", "phylogeny.jl"))
     include(joinpath(@__DIR__, "routes", "cladistics.jl"))
     include(joinpath(@__DIR__, "routes", "clade_cumulus.jl"))
-    include(joinpath(@__DIR__, "routes", "cladistics.jl"))
 
     ## Bad request middleware
     # Helpers deep in a route (e.g. `_req_group`) throw this to reject a request.

test/unit/test_clade_cumulus.jl

@@ -310,6 +310,10 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0
             con = DBInterface.connect(db)
             DuckDB.register_data_frame(con, df, "src_df")
             DBInterface.execute(con, "CREATE TABLE merged AS SELECT * FROM src_df")
+            # The registration is a view backed by a replacement scan that only
+            # exists on this connection, so it is removed again: left in the file
+            # it makes the view unreadable for every later connection.
+            DuckDB.unregister_data_frame(con, "src_df")
             DBInterface.close!(con)
             close(db)
             before_db_sha = bytes2hex(open(sha256, db_path))

test/unit/test_exact_summaries.jl

@@ -14,7 +14,10 @@ const EXACT_FEATURES = ["f1", "f2", "f3"]
     @test !s.approximate
     @test s.samples[1].total == 3
     @test s.samples[1].proportions == fill(big(1)//3, 3)
-    @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), nothing, nothing]
+    # Sample B has [2, 0, 0] reads: a zero count in a sample that has reads is a
+    # defined proportion of zero. Only a zero *total* leaves proportions
+    # undefined, and this sample's total is 2.
+    @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), big(0), big(0)]
     @test length(s.groups) == 2
     @test s.groups[1].label == "g1"
     @test s.groups[2].counts == BigInt[2, 0, 0]
@@ -28,7 +31,7 @@ const EXACT_FEATURES = ["f1", "f2", "f3"]
     @test storage["input_status"] == "exact"
     @test storage["samples"][1]["proportions"] == ["1/3", "1/3", "1/3"]
     @test storage["samples"][1]["total"] == 3
-    @test storage["samples"][2]["proportions"][2] === nothing
+    @test storage["samples"][2]["proportions"] == ["1/1", "0/1", "0/1"]
     @test length(jcs_sha256(storage)) == 64
     display = summary_to_display(s; digits=8)
     @test occursin("exact", display)

test/unit/test_validated_statistics.jl

@@ -92,16 +92,26 @@ end
 
 @testset "Numeric JSON representations fail closed" begin
     _vs_with_root() do _
-        for bad in (1.0, 1.5, -1, 9007199254740993)
+        # JSON3 normalises every integral numeric token - `1`, `1.0` and `1e0`
+        # alike - to Int64 while parsing, so those three are indistinguishable by
+        # the time a route sees them, and the endpoint must not refuse a token
+        # whose value is exactly the integer it denotes. What fails closed is
+        # everything that is not exactly an integer of Int64 size: a fractional
+        # token, a negative count and a token past the JSON-safe integer range.
+        for bad in (1.5, -1, 9007199254740993)
             payload = merge(_VS_EXACT_BODY, Dict("persist" => false,
                 "counts" => [[bad, 1], [1, 0], [1, 0]]))
             response = _vs_post("/api/v1/statistics/summaries", payload)
             @test response.status == 422
             @test JSON3.read(String(response.body)).error == "invalid_numeric_input"
         end
-        # Exponent notation is a float token even though its value is integral.
-        raw = """{"numeric_policy":{"mode":"exact_counts"},"sample_labels":["A"],"feature_labels":["f"],"counts":[[1e0]]}"""
-        @test _vs_post("/api/v1/statistics/summaries", raw; raw=true).status == 422
+        # 1e0 is an exactly integral token, so it is read as the count 1.
+        integral = """{"numeric_policy":{"mode":"exact_counts"},"sample_labels":["A"],"feature_labels":["f"],"counts":[[1e0]]}"""
+        @test _vs_post("/api/v1/statistics/summaries", integral; raw=true).status == 200
+        # 1.0000000001 is not integral, so the value that survives parsing is a
+        # Float64 and the route refuses it rather than rounding.
+        fractional = """{"numeric_policy":{"mode":"exact_counts"},"sample_labels":["A"],"feature_labels":["f"],"counts":[[1.0000000001]]}"""
+        @test _vs_post("/api/v1/statistics/summaries", fractional; raw=true).status == 422
         negative_string = merge(_VS_EXACT_BODY, Dict("persist" => false,
             "counts" => [["-1", 1], [1, 0], [1, 0]]))
         @test _vs_post("/api/v1/statistics/summaries", negative_string).status == 422