diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9bad3b7b..7b5a57a5 100644 --- a/.github/workflows/ci.yml +++ b/.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" diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index a702fde2..b4fa7247 100644 --- a/src/analysis/clade_cumulus.jl +++ b/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 diff --git a/src/pipeline/phylogeny.jl b/src/pipeline/phylogeny.jl index b394d9dc..b978ef02 100644 --- a/src/pipeline/phylogeny.jl +++ b/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 diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index 45825b61..eee7d490 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/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) diff --git a/src/server/routes/trees.jl b/src/server/routes/trees.jl index f3cc8e07..beec713a 100644 --- a/src/server/routes/trees.jl +++ b/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 diff --git a/src/server/routes/validated_statistics.jl b/src/server/routes/validated_statistics.jl index ce214c7d..082ef206 100644 --- a/src/server/routes/validated_statistics.jl +++ b/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 diff --git a/src/server/server.jl b/src/server/server.jl index ebf5abb6..eca368a0 100644 --- a/src/server/server.jl +++ b/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. diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index da09bc18..51d9d041 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/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)) diff --git a/test/unit/test_exact_summaries.jl b/test/unit/test_exact_summaries.jl index e1ace22d..d78f23e3 100644 --- a/test/unit/test_exact_summaries.jl +++ b/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) diff --git a/test/unit/test_validated_statistics.jl b/test/unit/test_validated_statistics.jl index c19bee2c..00a3214c 100644 --- a/test/unit/test_validated_statistics.jl +++ b/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