diff --git a/config/ci/lint_source.jl b/config/ci/lint_source.jl index fd11183c..ae50ba32 100644 --- a/config/ci/lint_source.jl +++ b/config/ci/lint_source.jl @@ -159,31 +159,118 @@ end # # `using Statistics` in Execution.jl with Statistics missing from [deps] made # precompilation fail outright. +# +# The first version of this check had two false negatives, and `import Printf` +# walked through both of them: it passed this gate locally and failed +# precompilation in CI with "Package MetaManifold does not have Printf in its +# dependencies" -- precisely the failure the check exists to prevent. +# +# 1. It exempted a list of "stdlibs that need no [deps] entry", Printf among +# them. There is no such class of name. A stdlib used by a package must be +# declared like any other dependency; only Base, Core and Main are bound +# without a declaration. Verified by experiment against every name the old +# list contained: each one fails with "does not have X in its dependencies" +# when undeclared. +# 2. It matched one name per line, so `using JSON3, YAML` checked JSON3 and +# never looked at YAML. +# +# Both are fixed by parsing the statement rather than pattern-matching its text: +# the AST has already resolved comma lists, relative imports and `using A: b, c` +# selection, none of which need guessing at. `self_test_declared_deps` below +# replays the exact input that defeated the old version, so a future rewrite that +# reintroduces either blind spot fails the gate instead of silently passing it. # --------------------------------------------------------------------------- -function check_declared_deps() - proj = read(joinpath(ROOT, "Project.toml"), String) - # `m` is required: without it `^` anchors to the whole string, not each line. - declared = Set(String[m.captures[1] for m in eachmatch(r"^([A-Za-z0-9_]+)\s*=\s*\""m, proj)]) + +# Bound in every module without a declaration. This list is deliberately tiny, +# and deliberately not a list of stdlibs: a stdlib is an ordinary dependency. +const ALWAYS_BOUND = Set(["Base", "Core", "Main"]) + +"""Root module of one `using`/`import` target, or `nothing` if it is relative.""" +function imported_root(x)::Union{Symbol,Nothing} + x isa Symbol && return x === :. ? nothing : x + x isa QuoteNode && return imported_root(x.value) + x isa Expr || return nothing + # `using A`, `using A.B`, `using A, B` -- first arg is the root, except for + # `using .A` / `using ..A`, which are relative and name no dependency. + x.head === :. && return isempty(x.args) ? nothing : imported_root(x.args[1]) + x.head === :(:) && return imported_root(x.args[1]) # `using A: b, c` + return nothing +end + +"""Package names `source` imports that `project_text` does not declare.""" +function undeclared_deps(project_text::AbstractString, source::AbstractString)::Vector{String} + declared = Set{String}(m.captures[1] for m in eachmatch(r"^([A-Za-z0-9_]+)\s*=\s*\""m, project_text)) # A package refers to itself by name inside its own source; that is not a dep. - self_name = match(r"^name\s*=\s*\"([A-Za-z0-9_]+)\""m, proj) + self_name = match(r"^name\s*=\s*\"([A-Za-z0-9_]+)\""m, project_text) self_name !== nothing && push!(declared, self_name.captures[1]) - # stdlibs that ship with Julia and need no [deps] entry - stdlib = Set(["Base", "Core", "Main", "Pkg", "Test", "UUIDs", "Dates", "Random", - "Printf", "Logging", "Statistics", "SHA", "Downloads", "LinearAlgebra", - "SparseArrays", "DelimitedFiles", "Sockets", "Markdown", "InteractiveUtils", - "Serialization", "Distributed", "Libdl", "Profile", "SuiteSparse"]) - for file in src_files() - for (i, line) in enumerate(eachline(file)) - s = strip(line) - startswith(s, '#') && continue - for m in eachmatch(r"^\s*(?:using|import)\s+([A-Za-z0-9_]+)(?![.\w])", s) - pkg = m.captures[1] - pkg in stdlib && continue - pkg in declared && continue - note("undeclared-dependency", file, i, - "`$pkg` is used but absent from Project.toml [deps]; precompilation will fail.") + + found = String[] + function walk(node) + node isa Expr || return nothing + if node.head === :using || node.head === :import + for arg in node.args + root = imported_root(arg) + root === nothing && continue + name = String(root) + (name in declared || name in ALWAYS_BOUND || name in found) && continue + push!(found, name) end end + for a in node.args + walk(a) + end + return nothing + end + walk(Meta.parseall(source)) + return found +end + +function check_declared_deps() + project = read(joinpath(ROOT, "Project.toml"), String) + for file in src_files() + lines = readlines(file) + for pkg in undeclared_deps(project, read(file, String)) + # Point at the statement, not at the top of the file. + ln = findfirst(l -> occursin(r"^\s*(?:using|import)\b", l) && occursin(pkg, l), lines) + note("undeclared-dependency", file, something(ln, 0), + "`$pkg` is used but absent from Project.toml [deps]; precompilation will fail.") + end + end +end + +# --------------------------------------------------------------------------- +# Check 4b — check 4, tested against the input that defeated it. +# +# A gate that cannot fail is decoration. Each case below is a defect that +# actually reached CI (the first two) or its inverse, which would be a false +# positive and just as damaging to trust in the gate. +# --------------------------------------------------------------------------- +function self_test_declared_deps() + fixture = """ + name = "Fixture" + + [deps] + JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + """ + cases = [ + # 1. The defect that motivated this rewrite: an undeclared STDLIB. + # `import Printf` passed the old check and failed CI precompilation. + ("using JSON3\nimport Printf\n", ["Printf"], "an undeclared stdlib"), + # 2. The second name in a comma list -- never examined by the old check. + ("using JSON3, YAML\n", ["YAML"], "a name after the first in a using list"), + # 3. Inverses: none of these is a dependency. + ("using ..Sibling\nusing .Local\nusing Fixture: Thing\nusing Fixture.Child\n", + String[], "relative and self imports"), + # 4. A declared stdlib is fine, and must not be reported. + ("import Statistics\n", String[], "a declared stdlib"), + ] + for (source, expected, what) in cases + got = sort(undeclared_deps(fixture, source)) + got == sort(expected) && continue + note("lint-self-test", joinpath("config", "ci", "lint_source.jl"), 0, + "the declared-dependency check is wrong about $what: expected $(sort(expected)), " * + "got $got. Fix the check before trusting it.") end end @@ -343,6 +430,7 @@ for (name, f) in [("parse", check_parses), ("escaped interpolation", check_escaped_interpolation), ("adjacent docstrings", check_adjacent_docstrings), ("declared dependencies", check_declared_deps), + ("dependency check self-test", self_test_declared_deps), ("test imports", check_test_imports)] n = length(failures) try diff --git a/docs/statistics/method-conditions/exact-descriptive-summaries.md b/docs/statistics/method-conditions/exact-descriptive-summaries.md index 1d15c6ed..66338b6c 100644 --- a/docs/statistics/method-conditions/exact-descriptive-summaries.md +++ b/docs/statistics/method-conditions/exact-descriptive-summaries.md @@ -8,6 +8,19 @@ SPDX-License-Identifier: CC-BY-SA-4.0 implementation is held to this document, so it is written to be checked against rather than admired. +**Status: implemented** in `src/analysis/exact_summaries.jl`. The conditions below were +**not** amended to fit the code; where the two disagreed, the code changed. Doing so +surfaced two defects in the layer underneath, which are fixed in the same change: +`to_display` printed a rendering labelled *6dp* with eighty digits after it, and it wrote +exact rationals as Julia's `2//3` — syntax leaking into text a person reads. + +**Evidence, as delivered** (in `test/unit/test_exact_summaries.jl`): hand-derived known +answers; an independent reference compared value by value against Python's +`fractions.Fraction` (skipping loudly by name if the runner has no `python3`); negative +controls for the zero-total case, for a float claiming exactness, for a negative count and +for a budget overrun; and a guard that the pipeline does not call this module, so +"unchanged when the layer is not selected" is a checked property rather than a promise. + ## What this method is Counts and proportions per sample — and per group when a grouping is supplied — carried at diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 1e86c791..c724aba7 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -29,6 +29,8 @@ include("pipeline/swarm.jl") # Analysis include("analysis/numeric_policy.jl") +# Exact summaries are catalogue item 1 and are built on the numeric policy, so they follow it. +include("analysis/exact_summaries.jl") include("analysis/diversity.jl") include("analysis/analysis.jl") include("analysis/AnalysisConfig.jl") diff --git a/src/analysis/exact_summaries.jl b/src/analysis/exact_summaries.jl new file mode 100644 index 00000000..3a964518 --- /dev/null +++ b/src/analysis/exact_summaries.jl @@ -0,0 +1,354 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Catalogue item 1 of docs/statistics/method-catalogue-v1.md: descriptive summaries at +# exact precision. +# +# The conditions this code is held to were published BEFORE it existed +# (docs/statistics/method-conditions/exact-descriptive-summaries.md). Where this file +# and that document disagree, the document is right and this file is the bug. +# +# The module exists to keep two claims apart that a single number cannot distinguish: +# +# a zero count is a value. A feature with no reads, in a sample that has reads, +# has a relative abundance of exactly 0. +# a zero total is not a small composition. A sample with no reads at all has NO +# relative abundances. Reporting 0 there does not round the truth, it +# invents a measurement -- and it is the kind of default that is only +# ever fixed before the code exists, which is why the document came +# first. +# +# No inference lives here. There is no p-value, no interval, no model and no comparison +# in this module's vocabulary -- not because they were left out, but because a +# descriptive summary that starts reporting significance is a different method wearing +# this one's name. +module ExactSummaries + +using ..NumericPolicy: NumericPolicySpec, numeric_policy, assert_mode, require_count, + checked_count_sum, exact_relative_abundance, exact_value, + to_storage, to_display, policy_fingerprint, + UnsupportedRepresentationError + +export ExactSummary, SampleSummary, GroupSummary, exact_summary, + has_defined_proportions, summary_to_storage, summary_to_display + +## Types + +""" + SampleSummary + +One sample's exact descriptive summary. + +`proportions[i]` is the exact relative abundance of `features[i]`, or `nothing` when +`total` is zero. The `nothing` is load-bearing: it is the difference between "this +feature was not seen" (exactly `0//1`) and "this sample has no composition to divide" +(undefined). They are not the same fact and no float can hold both. + +Counts are carried as integers of unbounded width. That is deliberate: the boundary +audit (#52) found real counts past 2^53, which no Float64 can hold, and normalising +them to a fixed width here would undo the audit at the last step. +""" +struct SampleSummary + label :: String + features :: Vector{String} + counts :: Vector{BigInt} + total :: BigInt + proportions :: Vector{Union{Rational{BigInt},Nothing}} + approximate :: Bool + warnings :: Vector{String} +end + +""" + GroupSummary + +Several samples aggregated exactly before being summarised. + +Aggregation happens in exact integer arithmetic, so a group total is a sum of counts +rather than a sum of proportions. It supports no comparison: a reader looking for a +difference between two of these should notice there is no place to put one. +""" +struct GroupSummary + label :: String + members :: Vector{String} + features :: Vector{String} + counts :: Vector{BigInt} + total :: BigInt + proportions :: Vector{Union{Rational{BigInt},Nothing}} + approximate :: Bool + warnings :: Vector{String} +end + +""" + ExactSummary + +The result of a descriptive summary: per-sample, per-group, and the policy it ran +under. + +The policy travels with the result because a summary that cannot say what numeric +policy produced it cannot be reproduced or compared with another run -- that is what +`policy_fingerprint` is for. +""" +struct ExactSummary + samples :: Vector{SampleSummary} + groups :: Vector{GroupSummary} + policy :: NumericPolicySpec + approximate :: Bool + warnings :: Vector{String} +end + +""" + has_defined_proportions(s::Union{SampleSummary,GroupSummary}) -> Bool + +False when the summary's total is zero, in which case its proportions are undefined +rather than zero. +""" +has_defined_proportions(s::Union{SampleSummary,GroupSummary}) = !iszero(s.total) + +## Construction + +# Reuse NumericPolicy's definition of what a count IS rather than restating it, and add +# the coordinates. "9007199254740994.0 is not a count" is a true sentence that gives a +# caller nothing to act on; naming the feature and sample is what makes it fixable. +function _exact_count(value, feature::AbstractString, sample::AbstractString) + count = try + require_count(value) + catch err + err isa UnsupportedRepresentationError || err isa ArgumentError || rethrow() + detail = err isa UnsupportedRepresentationError ? err.value : string(err) + throw(UnsupportedRepresentationError("an exactly-known count", string(value), + "feature '$feature', sample '$sample' ($detail)")) + end + count < 0 && throw(ArgumentError( + "negative count $count at feature '$feature', sample '$sample'; counts are " * + "non-negative, and a negative one here means the input is not a count table")) + return BigInt(count) +end + +function _summarise(label::String, features::Vector{String}, counts::Vector{BigInt}, + policy::NumericPolicySpec, approximate::Bool, + warnings::Vector{String}) + # `:widen` is defensive rather than load-bearing here: counts arrive from + # `_exact_count`, which already leaves them as BigInt, so this accumulator cannot + # overflow. It is kept because the alternative -- a future refactor that carries + # Int64 for speed -- would otherwise wrap silently, and a wrapped total is a wrong + # number no downstream check can detect. Mutation testing confirmed the flag is not + # exercised today; the test that matters asserts the total, not the mechanism. + total = checked_count_sum(counts; on_overflow = :widen) + proportions = Vector{Union{Rational{BigInt},Nothing}}(undef, length(counts)) + if iszero(total) + # every feature's proportion is undefined, and each one says so by being nothing + fill!(proportions, nothing) + push!(warnings, "'$label' has a total of zero reads: its relative abundances " * + "are undefined and are carried as nothing, never as 0 -- a " * + "proportion that does not exist is not a proportion of zero") + else + for (i, count) in enumerate(counts) + proportions[i] = exact_relative_abundance(count, total; policy) + end + end + return (total, proportions, warnings) +end + +function _sample_summary(label::String, features::Vector{String}, + counts::Vector{BigInt}, policy::NumericPolicySpec, + approximate::Bool) + warnings = String[] + approximate && push!(warnings, + "'$label' was supplied as floating point: its values are carried as " * + "approximations and labelled as such, because a Float64 above 2^53 cannot say " * + "which integer it holds and that history cannot be inspected") + total, proportions, warnings = _summarise(label, features, counts, policy, + approximate, warnings) + return SampleSummary(label, features, counts, total, proportions, approximate, warnings) +end + +function _group_summary(label::String, members::Vector{String}, features::Vector{String}, + row_counts::Vector{Vector{BigInt}}, policy::NumericPolicySpec, + approximate::Bool) + warnings = String[] + approximate && push!(warnings, + "'$label' aggregates samples supplied as floating point: its values are " * + "carried as approximations and labelled as such") + # Sums, not averaged proportions: a proportion averaged across unequal depths + # weights a shallow sample the same as a deep one, which is a different (and + # usually unintended) statement. + counts = BigInt[checked_count_sum(row; on_overflow = :widen) for row in row_counts] + total, proportions, warnings = _summarise(label, features, counts, policy, + approximate, warnings) + return GroupSummary(label, members, features, counts, total, proportions, + approximate, warnings) +end + +""" + exact_summary(counts; sample_labels, feature_labels, groups, policy) -> ExactSummary + +Exactly-summarise a features-by-samples count table. + +Rows are features and columns are samples, matching the rest of the analysis layer. The +table's values must be counts: integers, or floats that are finite, integral and within +2^53 (a float beyond that is refused, because it cannot say which integer it holds). +Float input is accepted but the result is marked `approximate` and says so in its +warnings -- approximate input is not hidden, it is labelled. + +`groups`, when given, names one group per sample; the group summary aggregates the +member samples in exact integer arithmetic first. + +The policy must be `:exact_counts`. `:ordinary` is refused by name, because a caller +who asks for an exact summary under an ordinary policy would otherwise receive +Float64s that look exactly like the exact answer until someone checks; and +`:high_precision` is refused too, since higher precision is not exactness -- it moves +the rounding error rather than abolishing it (see `numeric-contracts.md`). + + julia> s = exact_summary([4 6; 0 3]; + sample_labels = ["a", "b"], + feature_labels = ["f1", "f2"]); + + julia> s.samples[2].proportions + 2-element Vector{Union{Nothing, Rational{BigInt}}}: + 2//3 + 1//3 +""" +function exact_summary(counts::AbstractMatrix{<:Real}; + sample_labels::AbstractVector{<:AbstractString}, + feature_labels::AbstractVector{<:AbstractString} = + ["feature $i" for i in 1:size(counts, 1)], + groups::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + policy::NumericPolicySpec = numeric_policy(:exact_counts)) + assert_mode(policy, :exact_counts) + + n_features, n_samples = size(counts) + length(sample_labels) == n_samples || throw(ArgumentError( + "expected $n_samples sample labels for $n_samples columns, got " * + "$(length(sample_labels))")) + length(feature_labels) == n_features || throw(ArgumentError( + "expected $n_features feature labels for $n_features rows, got " * + "$(length(feature_labels))")) + + features = String[String(f) for f in feature_labels] + labels = String[String(s) for s in sample_labels] + + approximate = !(eltype(counts) <: Integer) + warnings = String[] + + columns = Vector{Vector{BigInt}}(undef, n_samples) + for j in 1:n_samples + columns[j] = BigInt[_exact_count(counts[i, j], features[i], labels[j]) + for i in 1:n_features] + end + + samples = [_sample_summary(labels[j], features, columns[j], policy, approximate) + for j in 1:n_samples] + + group_summaries = GroupSummary[] + if !isnothing(groups) + length(groups) == n_samples || throw(ArgumentError( + "expected $n_samples group labels (one per sample), got $(length(groups))")) + group_labels = String[String(g) for g in groups] + for label in unique(group_labels) + members = [labels[j] for j in 1:n_samples if group_labels[j] == label] + row_counts = [BigInt[columns[j][i] for j in 1:n_samples + if group_labels[j] == label] for i in 1:n_features] + push!(group_summaries, _group_summary(String(label), members, features, + row_counts, policy, approximate)) + end + end + + return ExactSummary(samples, group_summaries, policy, approximate, warnings) +end + +## Boundaries + +# Storage and display are separated deliberately: what is written down must round-trip +# exactly, and what a person reads may be rendered. The rendered form is never the +# stored one, which is the whole reason `to_storage` exists on the numeric layer. +function _counts_to_storage(x::Union{SampleSummary,GroupSummary}) + return Dict{String,Any}( + "label" => x.label, + "total" => to_storage(exact_value(x.total)), + "counts" => [to_storage(exact_value(c)) for c in x.counts], + "proportions" => [isnothing(p) ? nothing : to_storage(exact_value(p)) + for p in x.proportions], + "proportions_defined" => has_defined_proportions(x), + "approximate" => x.approximate, + "warnings" => x.warnings, + ) +end + +""" + summary_to_storage(s::ExactSummary) -> Dict{String,Any} + +The summary as it should be written down: exact counts and proportions as integers or +as `"numerator/denominator"` strings, `nothing` where a proportion is undefined, and +the policy fingerprint so the run can be compared with another. + +A proportion stored here parses back with `parse_exact_rational` to the identical +rational; a rounded decimal never appears in this structure. +""" +function summary_to_storage(s::ExactSummary) + out = Dict{String,Any}( + "kind" => "descriptive summary", + "claim" => "counts and proportions are exact; no comparison, test or " * + "significance is claimed", + "policy_fingerprint" => policy_fingerprint(s.policy), + "mode" => String(s.policy.mode), + "approximate" => s.approximate, + "samples" => [_counts_to_storage(x) for x in s.samples], + "groups" => [_counts_to_storage(x) for x in s.groups], + "warnings" => s.warnings, + ) + for (i, x) in enumerate(s.samples) + out["samples"][i]["features"] = x.features + end + for (i, x) in enumerate(s.groups) + out["groups"][i]["features"] = x.features + out["groups"][i]["members"] = x.members + end + return out +end + +""" + summary_to_display(s::ExactSummary; digits=6) -> String + +The summary as a person reads it: exact fractions, with a rounded decimal offered +beside them and marked as a rendering. + +The first line states what the numbers are not. A table of per-group counts is the +shape most often misread as a comparison, so the text says once, plainly, that no +comparison was made. +""" +function summary_to_display(s::ExactSummary; digits::Integer = 6) + digits >= 0 || throw(ArgumentError("digits must not be negative, got $digits")) + io = IOBuffer() + println(io, "Descriptive summary — counts and proportions are exact. No comparison, " * + "no test and no significance is claimed.") + if s.approximate + println(io, " ⚠ some input was floating point: those values are approximations, " * + "not exact counts") + end + for (heading, entries) in (("sample", s.samples), ("group", s.groups)) + for x in entries + println(io, " $heading ", x.label, " — total ", + to_display(exact_value(x.total); digits = digits)) + if !has_defined_proportions(x) + println(io, " proportions: undefined (zero total) — not zero") + continue + end + for (feature, proportion) in zip(x.features, x.proportions) + println(io, " ", feature, ": ", + to_display(exact_value(proportion); digits = digits)) + end + end + end + # Built by hand rather than vcat'ed with a generator: `vcat` treats a generator as a + # scalar here, and the note line printed the iterator's own type instead of a warning. + notes = copy(s.warnings) + for entry in vcat(s.samples, s.groups), warning in entry.warnings + push!(notes, warning) + end + for warning in notes + println(io, " note: ", warning) + end + return String(take!(io)) +end + +end # module ExactSummaries diff --git a/src/analysis/numeric_policy.jl b/src/analysis/numeric_policy.jl index 5798a378..2974c77c 100644 --- a/src/analysis/numeric_policy.jl +++ b/src/analysis/numeric_policy.jl @@ -433,6 +433,34 @@ function to_storage(v::NumericValue) return x isa BigFloat ? string(x) : Float64(x) end +# The rendering at `digits` places, computed in INTEGER arithmetic. +# +# The first version used BigFloat and a format string, which fixes the symptom -- a +# label saying 6dp followed by eighty digits, because `round(x; digits)` keeps the +# significand -- but adds a dependency to do it, and float arithmetic to render a value +# whose whole point is that it is not a float. This version divides in integers: scaling +# by 10^digits and rounding half away from zero cannot lose the value, works for numbers +# far beyond Float64, and needs nothing outside Base. +function _rounded_scaled(num::Integer, den::Integer, scale::Integer) + n = big(num) * scale + q, r = divrem(n, den) + if 2 * abs(r) >= abs(den) + q += sign(n) * sign(den) # half away from zero, in both directions + end + return q +end + +function _rendered_decimal(x::Rational, digits::Integer) + sign_str = numerator(x) < 0 ? "-" : "" + if digits == 0 + return string(sign_str, _rounded_scaled(abs(numerator(x)), denominator(x), big(1))) + end + scale = big(10)^digits + scaled = _rounded_scaled(abs(numerator(x)), denominator(x), scale) + whole, frac = divrem(abs(scaled), scale) + return string(sign_str, whole, ".", lpad(string(frac), digits, "0")) +end + """ to_display(v::NumericValue; digits=6) -> String @@ -444,10 +472,14 @@ before it arrives. function to_display(v::NumericValue; digits::Integer = 6) digits >= 0 || throw(ArgumentError("digits must not be negative, got $digits")) if v.kind === EXACT + # An exact rational renders as "2/3", not Julia's "2//3": the double slash is + # syntax leaking into something a person reads, and it also means the fraction + # shown is not the fraction that may be copied out. to_storage already writes + # numerator/denominator, so display and storage now agree on the form. return v.value isa Integer ? string(v.value) : - string(v.value, " (exact; ", digits, "dp = ", - round(BigFloat(numerator(v.value)) / BigFloat(denominator(v.value)); - digits = digits), ")") + string(numerator(v.value), "/", denominator(v.value), + " (exact; ", digits, "dp = ", + _rendered_decimal(v.value, digits), ")") end kept = v.kind === ROUNDED && !isnothing(v.digits) ? v.digits : Int(digits) return string(round(v.value; digits = kept), " (", v.kind, ", ", kept, "dp)") diff --git a/test/runtests.jl b/test/runtests.jl index 181b2a4c..5305f1f4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,6 +35,7 @@ using MetaManifold: AnalysisConfig include("unit/test_numeric_policy.jl") include("unit/test_numeric_boundaries.jl") + include("unit/test_exact_summaries.jl") include("unit/test_diversity.jl") include("unit/test_merge_taxa.jl") include("unit/test_config.jl") diff --git a/test/unit/test_exact_summaries.jl b/test/unit/test_exact_summaries.jl new file mode 100644 index 00000000..34d39bf2 --- /dev/null +++ b/test/unit/test_exact_summaries.jl @@ -0,0 +1,381 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# Evidence for catalogue item 1, held to the conditions published in +# docs/statistics/method-conditions/exact-descriptive-summaries.md before the +# implementation existed. Each testset below is one of the four things that document +# requires: known answers derived by hand, an independent reference, negative controls, +# and proof that nothing outside this module changed. + +using Test +using JSON3 +using MetaManifold.NumericPolicy +using MetaManifold.ExactSummaries + +const REPO_ROOT = normpath(joinpath(@__DIR__, "..", "..")) +const EXECUTION_SOURCE = joinpath(REPO_ROOT, "src", "analysis", "Execution.jl") + +@testset "exact descriptive summaries" begin + + @testset "known answers, derived by hand" begin + # features x samples: f1 = (4, 6), f2 = (0, 3) + counts = [4 6; + 0 3] + summary = exact_summary(counts; + sample_labels = ["a", "b"], + feature_labels = ["f1", "f2"]) + + @test summary.approximate == false + @test [s.label for s in summary.samples] == ["a", "b"] + + # Worked by hand, not read off the output: + # a: total 4 -> f1 = 4/4 = 1, f2 = 0/4 = 0 + # b: total 9 -> f1 = 6/9 = 2/3, f2 = 3/9 = 1/3 + a, b = summary.samples + @test a.total == 4 + @test a.proportions == [big(1) // big(1), big(0) // big(1)] + @test b.total == 9 + @test b.proportions == [big(2) // big(3), big(1) // big(3)] + + # A zero count in a sample that has reads is exactly zero -- a value, not a + # missing one. The distinction from the zero-total case below is the point. + @test !isnothing(a.proportions[2]) + @test a.proportions[2] == 0 # exactly zero, and provably a Rational + @test has_defined_proportions(a) + + # The representation itself is the claim: a Float64 here would round, and the + # rounding would be invisible downstream. + @test all(p -> p isa Rational{BigInt}, b.proportions) + + # Thirds and fifths have no exact binary form, so the exact values cannot be + # reproduced by the float that would normally carry them. + thirds = exact_summary(reshape([1, 2, 2], 3, 1); + sample_labels = ["c"], + feature_labels = ["g1", "g2", "g3"]).samples[1] + @test thirds.total == 5 + @test thirds.proportions == [big(1) // big(5), big(2) // big(5), big(2) // big(5)] + @test sum(thirds.proportions) == big(1) // big(1) # exact, not 0.9999... + @test Float64(big(1) // big(3)) != big(1) // big(3) + end + + @testset "counts past 2^53 survive exactly" begin + # The probe is the smallest integer a Float64 cannot hold, from the boundary + # audit (#52). Here it is a count, and the summary is the reason the audit + # mattered: this number cannot pass through a float, so it must not have to. + huge = big(2)^53 + 1 + summary = exact_summary(reshape([huge, big(1)], 2, 1); + sample_labels = ["deep"], + feature_labels = ["f1", "f2"]) + deep = summary.samples[1] + + @test deep.total == huge + 1 + @test deep.proportions[1] == huge // (huge + 1) + @test deep.proportions[2] == big(1) // (huge + 1) + @test numerator(deep.proportions[1]) == huge + + # The float claim fails, which is the whole reason for the string boundary. + @test Float64(huge) == big(2)^53 + + # A total beyond Int64 is widened, not wrapped: a wrapped total is a wrong + # number that no downstream check can detect. + # + # The inputs here are Int64 ON PURPOSE, and the trap is real: Julia's own sum + # wraps on exactly this fixture, which the assertion below demonstrates rather + # than asserts in prose. An implementation that adds these counts naively returns + # a negative total, and the assertion above fails. + # + # What this test does NOT claim: that the widening flag in the module is what + # saves it. Counts leave the boundary as BigInt, so the accumulator cannot + # overflow either way -- removing the flag still passes, which mutation testing + # showed. The claim under test is the total, not the mechanism; a mechanism test + # here would be testing an implementation choice. + @test sum([typemax(Int64), 1]) < 0 # the naive sum wraps: the trap is real + @test sum([typemax(Int64), 1]) != big(2)^63 + overflow = exact_summary(reshape([typemax(Int64), 1], 2, 1); + sample_labels = ["wide"], feature_labels = ["f1", "f2"]) + @test overflow.samples[1].total == big(2)^63 + @test overflow.samples[1].total > 0 + @test overflow.samples[1].proportions == + [(big(2)^63 - 1) // big(2)^63, big(1) // big(2)^63] + end + + @testset "a zero-total sample has no proportions (negative control)" begin + counts = [0 5; + 0 0] + summary = exact_summary(counts; + sample_labels = ["empty", "fine"], + feature_labels = ["f1", "f2"]) + empty, fine = summary.samples + + @test empty.total == 0 + @test has_defined_proportions(empty) == false + @test all(isnothing, empty.proportions) + + # The control: not zero, not NaN, not a float. A proportion that does not exist + # must not be representable as a number at all. + @test !any(p -> p == 0, empty.proportions) + @test !any(p -> p isa AbstractFloat, empty.proportions) + + # And it says which sample, in words, without anyone having to infer it. + @test any(w -> occursin("empty", w) && occursin("zero", lowercase(w)), empty.warnings) + + # A sample with reads in the same table is unaffected: the undefined state is + # per sample, not a property of the table. + @test has_defined_proportions(fine) + @test fine.proportions == [big(1) // big(1), big(0) // big(1)] + + # An all-zero table is a real observation, not an error. + allzero = exact_summary(zeros(Int, 2, 1); + sample_labels = ["z"], feature_labels = ["f1", "f2"]) + @test allzero.samples[1].total == 0 + @test all(isnothing, allzero.samples[1].proportions) + end + + @testset "per-group summaries aggregate exactly, and compare nothing" begin + counts = [4 6 1; + 0 3 2] + summary = exact_summary(counts; + sample_labels = ["a", "b", "c"], + feature_labels = ["f1", "f2"], + groups = ["A", "B", "A"]) + + @test length(summary.groups) == 2 + A = summary.groups[1] + @test A.label == "A" + @test A.members == ["a", "c"] # exact members, named + @test A.counts == [BigInt(5), BigInt(2)] # 4+1 and 0+2, in integers + @test A.total == 7 + @test A.proportions == [big(5) // big(7), big(2) // big(7)] + + # Aggregation is a sum of counts, not a mean of proportions: with depths 4 and 3 + # the two differ, and only one of them is the group's composition. + mean_of_proportions = (big(1) // big(1) + big(1) // big(3)) / 2 + @test A.proportions[1] != mean_of_proportions + + # A group whose members have no reads at all is undefined, exactly like a + # zero-total sample -- grouping must not quietly manufacture a composition. + zero_group = exact_summary(counts; + sample_labels = ["a", "b", "c"], + feature_labels = ["f1", "f2"], + groups = ["A", "B", "A"]).groups[1] + @test has_defined_proportions(zero_group) + + empty_group = exact_summary([0 1; 0 2]; + sample_labels = ["x", "y"], + feature_labels = ["f1", "f2"], + groups = ["E", "F"]).groups[1] + @test empty_group.total == 0 + @test !has_defined_proportions(empty_group) + @test all(isnothing, empty_group.proportions) + + # The shape people misread as a test says, in words, that it is not one. + text = summary_to_display(summary) + @test occursin("No comparison", text) + @test occursin("no significance is claimed", text) + end + + @testset "input that is not a count is refused by name (negative control)" begin + # A non-integral value is not a count, and the refusal names the cell. + err = try + exact_summary(reshape([1.5, 2.0], 2, 1); + sample_labels = ["s1"], feature_labels = ["f1", "f2"]) + nothing + catch e + e + end + @test err isa UnsupportedRepresentationError + @test occursin("f1", sprint(showerror, err)) + @test occursin("s1", sprint(showerror, err)) + + # A float beyond 2^53 cannot say which integer it holds -- even when it is a + # value that exists (2^53 + 2 is representable). Its history cannot be + # inspected, so it is refused rather than trusted. + @test_throws UnsupportedRepresentationError exact_summary( + reshape([2.0^53 + 2, 1.0], 2, 1); + sample_labels = ["s1"], feature_labels = ["f1", "f2"]) + + # Negative counts are not counts of a table; that is a broken input, not a + # rounding question. + @test_throws ArgumentError exact_summary([-1 2]; + sample_labels = ["s1"], + feature_labels = ["f1"]) + + # Float input inside the exact range is ACCEPTED, and labelled: the conditions + # allow an approximation to travel, never to travel unmarked. + approximate = exact_summary([4.0 6.0; + 0.0 3.0]; + sample_labels = ["a", "b"], + feature_labels = ["f1", "f2"]) + @test approximate.approximate == true + @test any(w -> occursin("approximation", w), approximate.samples[1].warnings) + @test approximate.samples[2].proportions == [big(2) // big(3), big(1) // big(3)] + @test occursin("approximations", summary_to_display(approximate)) + end + + @testset "the budget raises instead of rounding (negative control)" begin + # A prime total whose denominator cannot fit a 32-bit budget: the summary must + # fail loudly rather than fall back to floats, because a silent downgrade from + # exact to approximate is the failure this whole layer exists to prevent. + tight = numeric_policy(; mode = :exact_counts, max_denominator_bits = 32) + huge_prime = big(2)^61 - 1 + @test_throws ResourceLimitError exact_summary( + reshape([huge_prime, big(1)], 2, 1); + sample_labels = ["s"], feature_labels = ["f1", "f2"], policy = tight) + + # The same table under the default budget is fine: the limit is a budget, not a + # defect in the data. + generous = exact_summary(reshape([huge_prime, big(1)], 2, 1); + sample_labels = ["s"], feature_labels = ["f1", "f2"]) + @test generous.samples[1].proportions[1] == huge_prime // (huge_prime + 1) + end + + @testset "an ordinary policy cannot produce an exact summary" begin + # Asking for exactness under the ordinary policy must fail by name. The failure + # this prevents is a caller receiving Float64s that look exactly like the exact + # answer would have looked, until somebody checks. + err = try + exact_summary([4 6; 0 3]; + sample_labels = ["a", "b"], feature_labels = ["f1", "f2"], + policy = numeric_policy(:ordinary)) + nothing + catch e + e + end + @test err isa UnsupportedRepresentationError + @test occursin("exact_counts", sprint(showerror, err)) + + # Higher precision is not exactness: it moves the rounding error rather than + # abolishing it, so it cannot be accepted here either. + @test_throws UnsupportedRepresentationError exact_summary( + [4 6; 0 3]; sample_labels = ["a", "b"], feature_labels = ["f1", "f2"], + policy = numeric_policy(; mode = :high_precision, precision_bits = 256)) + end + + @testset "display renders signs and ties by integer arithmetic" begin + # The rendering exists because a rounded decimal must never be what gets stored, + # but it is still text a person reads, so it has to be right at the edges. The + # implementation does this in integers; these are the cases a floating-point + # renderer gets wrong, which is why they are here rather than assumed. + @testset "negative rationals keep their sign, and only one" begin + text = to_display(exact_value(-2 // 3)) + @test occursin("-0.666667", text) + @test !occursin("--", text) + end + @testset "a true tie rounds away from zero in both directions" begin + # 1/8 = 0.125 exactly, so at two decimals the tie is real and visible. + @test occursin("0.13", to_display(exact_value(1 // 8); digits = 2)) + @test occursin("-0.13", to_display(exact_value(-1 // 8); digits = 2)) + end + @testset "values beyond Float64 still render as digits" begin + # Float64 overflows its exact-integer range around 2^53; this is far past it, + # so a float-based renderer would lose the integer part or switch to e+00. + huge = (big(2)^200 + 1) // big(3) + text = to_display(exact_value(huge)) + @test occursin(".666667", text) + @test !occursin("e+", text) + @test !occursin("E+", text) + end + end + + @testset "storage round-trips through the boundary readers" begin + summary = exact_summary([4 6; 0 3]; + sample_labels = ["a", "b"], feature_labels = ["f1", "f2"]) + stored = summary_to_storage(summary) + + @test stored["claim"] isa String + @test occursin("no comparison", stored["claim"]) + @test stored["mode"] == "exact_counts" + @test length(stored["policy_fingerprint"]) == 64 + + b = stored["samples"][2] + @test b["total"] == 9 + # Proportions travel as strings, never as JSON numbers: a JSON number on the way + # through a browser is a Float64 (see the boundary audit). + @test b["proportions"] == ["2/3", "1/3"] + @test b["proportions_defined"] == true + + # The round trip: what was written parses back to the identical rational, so the + # stored form loses nothing. + @test NumericPolicy.parse_exact_rational(b["proportions"][1]) == big(2) // big(3) + @test NumericPolicy.parse_exact_rational(b["proportions"][2]) == big(1) // big(3) + + # Undefined travels as null, and says it is undefined rather than zero. + empty = summary_to_storage(exact_summary([0 5; 0 0]; + sample_labels = ["empty", "fine"], + feature_labels = ["f1", "f2"])) + @test empty["samples"][1]["proportions"] == [nothing, nothing] + @test empty["samples"][1]["proportions_defined"] == false + + # Display and storage are different things: the rendering is marked as one. + text = summary_to_display(summary) + @test occursin("2/3 (exact; 6dp = 0.666667)", text) # exact fraction, rendering beside it + @test !occursin("0.6666669999", text) # the rendering is rounded, not 80 digits + @test !occursin("2/3 (exact;", b["proportions"][1]) # never what gets stored + end + + @testset "independent reference: python3 fractions" begin + # Compared value by value against an independent implementation of exact rational + # arithmetic. Python's fractions.Fraction reduces to lowest terms by its own code, + # so agreement is evidence about the arithmetic rather than about this module's + # agreement with itself. + python = Sys.which("python3") + if isnothing(python) + @testset "python3 absent — the independent reference did NOT run" begin + @test_skip false + end + else + counts = [7 3; + 5 11; + 13 2] + summary = exact_summary(counts; + sample_labels = ["a", "b"], + feature_labels = ["f1", "f2", "f3"]) + # Joined from lines rather than written as an indented triple-quoted + # string: a stray leading space on the first line is an IndentationError in + # Python, and the failure then looks like a broken reference rather than a + # typo in the test. + script = join([ + "import json", + "from fractions import Fraction", + "counts = [[7, 3], [5, 11], [13, 2]]", + "out = []", + "for j in range(len(counts[0])):", + " col = [counts[i][j] for i in range(len(counts))]", + " total = sum(col)", + " out.append({\"total\": total, \"proportions\": " * + "[str(Fraction(c, total).numerator) + \"/\" + " * + "str(Fraction(c, total).denominator) for c in col]})", + "print(json.dumps(out))", + ], "\n") + reference = JSON3.read(read(`$python -c $script`, String)) + + for (j, sample) in enumerate(summary.samples) + @test sample.total == reference[j]["total"] + for (i, proportion) in enumerate(sample.proportions) + mine = string(numerator(proportion), "/", denominator(proportion)) + @test mine == reference[j]["proportions"][i] + end + end + for sample in summary.samples + @test sum(sample.proportions) == big(1) // big(1) + end + end + end + + @testset "the existing pipeline does not depend on this module" begin + # "Unchanged behaviour when the layer is not selected" is satisfied by the layer + # genuinely being additional: nothing in the default path calls it. Wiring it into + # the pipeline would change results that saved analyses were computed from, which + # is a decision, so it fails here by name instead of happening quietly. + execution = read(EXECUTION_SOURCE, String) + @test !occursin("ExactSummaries", execution) + @test !occursin("exact_summary", execution) + + # Running a summary does not touch the ambient policy: the ordinary default is + # still the default afterwards, which is what keeps concurrent analyses from + # sharing a setting. + before = numeric_policy() + @test before.mode === :ordinary + exact_summary([4 6; 0 3]; sample_labels = ["a", "b"], feature_labels = ["f1", "f2"]) + @test numeric_policy().mode === :ordinary + @test policy_fingerprint(numeric_policy()) == policy_fingerprint(numeric_policy(:ordinary)) + end +end