Resolve PR #224 against main: validated-statistics + cladistics fixes, review findings, and merge conflicts

.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"
 

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

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/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

@@ -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

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

@@ -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))

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)

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