diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d51aa441..5fa55493 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -400,11 +400,26 @@ 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 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 failed (tail)::$(tail -c 3500 /tmp/julia-validated-tests.log | enc)" + enc() { head -c 3000 | 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 2500 /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 -c 3000 /tmp/validated-statistics-tests.log | enc)" fi exit "$status" diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 63cbdf24..9ac3f4f0 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -28,11 +28,15 @@ const DEFAULT_TAXONOMIC_RANKS = [ "Class", "Order", "Family", "Genus", "Species", ] -# Placement tools may truncate a query's candidate list, so retained -# like-weight ratios can legitimately sum to slightly less than 1 (a single -# retained placement at 0.98 is complete evidence). Sums far from 1 indicate -# a malformed or partial record and stay undecided. +# Placement tools truncate a query's candidate list at a likelihood-weight +# cutoff, so retained like-weight ratios can legitimately sum to a little +# under 1 (a single retained placement at 0.98 is complete evidence). A shortfall +# beyond the tolerance is a missing or discarded placement and stays undecided. +# A sum over 1 is never legitimate: weights are mutually exclusive placement +# ratios, so overweight evidence is malformed and is rejected by a tight upper +# bound while the lower shortfall is tolerated. const PLACEMENT_WEIGHT_SUM_TOL = 0.05 +const PLACEMENT_WEIGHT_OVER_TOL = 1e-9 ## Phylogenetic tree structures & Newick / jplace parsers mutable struct PhyloNode @@ -573,7 +577,8 @@ function parse_jplace(raw_json::AbstractString)::JplaceDocument if !isempty(bad_reason) push!(queries, QueryPlacement(q_ord, q_names, QueryCandidate[], :undecided, bad_reason, nothing)) - elseif abs(weight_sum - 1.0) > PLACEMENT_WEIGHT_SUM_TOL + elseif weight_sum > 1.0 + PLACEMENT_WEIGHT_OVER_TOL || + weight_sum < 1.0 - PLACEMENT_WEIGHT_SUM_TOL push!(queries, QueryPlacement(q_ord, q_names, cands, :undecided, "weight_sum_out_of_tolerance ($(round(weight_sum; digits=6)))", nothing)) else 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/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 fb14d584..e6cf4ffd 100644 --- a/src/server/routes/validated_statistics.jl +++ b/src/server/routes/validated_statistics.jl @@ -310,8 +310,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 @@ -322,6 +320,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[]) @@ -348,45 +355,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], ", ") - # Rows sharing a feature (e.g. several sequences in one genus) are summed - # here, in BigInt, not in SQL: DuckDB refuses bare count columns beside - # GROUP BY, and SUM over HUGEINT counts would overflow exactly the values - # this path exists to keep exact. ORDER BY keeps features sorted and each - # feature's rows adjacent. + # 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 ORDER BY feature LIMIT $(MAX_CELLS + 1)" frame = DataFrame(DBInterface.execute(con, sql, where_params)) - nrow(frame) <= MAX_CELLS || throw(ResourceLimitError("rows", MAX_CELLS, nrow(frame))) nrow(frame) > 0 || throw(ArgumentError("no rows remain after filtering")) + nrow(frame) <= MAX_CELLS || throw(ResourceLimitError("table rows", MAX_CELLS, nrow(frame))) features = String[] - feature_rows = Vector{BigInt}[] + 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")) - feature = string(row.feature) - if isempty(features) || features[end] != feature - push!(features, feature) - push!(feature_rows, zeros(BigInt, length(samples))) + 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 - target = feature_rows[end] for (j, sample) in enumerate(samples) value = row[Symbol(sample)] ismissing(value) && throw(UnsupportedRepresentationError( - "a non-missing integer count", "NULL", "feature '$feature', 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 '$feature', sample '$sample'")) + "an integer count", repr(value), "feature '$label', sample '$sample'")) value >= 0 || throw(ArgumentError( - "negative count $value at feature '$feature', sample '$sample'")) - target[j] += BigInt(value) # every addition 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, feature_row) in enumerate(feature_rows) - counts[i, :] = feature_row + for (i, label) in enumerate(features) + counts[i, :] = totals[label] end - length(features) * length(samples) <= MAX_CELLS || throw(ResourceLimitError( - "table cells", MAX_CELLS, length(features) * length(samples))) 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 6e12a0c9..51d9d041 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -252,8 +252,8 @@ done echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_000004:0.08);" > "$OUT" """) chmod(mock_ft, 0o755) - old_ft = get(MetaManifold.Tools._tools, "fasttree", nothing) - MetaManifold.Tools._tools["fasttree"] = mock_ft + old_ft = get(PH.Tools._tools, "fasttree", nothing) + PH.Tools._tools["fasttree"] = mock_ft try s4 = [ "ASV_1" => "ACGTACGTACGTACGTACGT", @@ -267,9 +267,9 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 r["sha256"] == bytes2hex(open(sha256, mock_ft)), ft_res.tool_records) finally if old_ft === nothing - delete!(MetaManifold.Tools._tools, "fasttree") + delete!(PH.Tools._tools, "fasttree") else - MetaManifold.Tools._tools["fasttree"] = old_ft + PH.Tools._tools["fasttree"] = old_ft end end end @@ -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 098667c8..d78f23e3 100644 --- a/test/unit/test_exact_summaries.jl +++ b/test/unit/test_exact_summaries.jl @@ -14,10 +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) - # A zero count in a sample with a positive total has an exact proportion - # of 0 (the independent Python fractions reference below agrees); only a - # zero total leaves proportions undefined. - @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1)//1, big(0)//1, big(0)//1] + # 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] @@ -31,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] == "0/1" + @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 7a7465d2..bea38a57 100644 --- a/test/unit/test_validated_statistics.jl +++ b/test/unit/test_validated_statistics.jl @@ -106,6 +106,10 @@ end raw = """{"numeric_policy":{"mode":"exact_counts"},"sample_labels":["A"],"feature_labels":["f"],"counts":[[$token]]}""" @test _vs_post("/api/v1/statistics/summaries", raw; raw=true).status == 422 end + + # Fractional input stays invalid as well, even when it is nearly 1. + 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