From 4e07b5a02e8ff97c03e0f13dfab9419194360b01 Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Mon, 25 May 2026 18:58:13 +0200 Subject: [PATCH 01/11] Fix correctness bugs surfaced in pre-JOSS code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Addresses ten severe findings from the pre-submission review: - StaticDesigns: fix UndefVarError (`nrow(data)` → `nrow(X)`) in the no-target branch. - fronts: rewrite `front()` with a proper Pareto-dominance predicate over joint minimization and a separate tolerance-thinning pass; previous implementation kept dominated points and dropped non-dominated ones. - Mahalanobis: canonicalize evidence-key ordering against the column order used to precompute Λ, so distances no longer depend on evidence insertion order. - EfficientValueMDP: tighten `value` `hasmethod` check to `Tuple{Evidence, NTuple{2, Float64}}` to match the documented signature. - efficient_designs / conditional_efficient_designs: guard `thresholds < 2` with an `ArgumentError` before calling `range(0.0, 1.0, thresholds)`. - POMDPs.actions: drop a multi-feature experiment iff ANY of its features is already observed (atomic semantics) across all three MDPs. - Variance() / Entropy(): throw `ArgumentError` when the prior gives zero initial variance/entropy instead of producing NaN ratios. - DistanceBased weights: extract `apply_masks!`; the zero-similarity recovery branch now re-applies hard-match and importance/filter masks, and throws if no historical row satisfies the constraints. - Public design entry-points (`efficient_design`, `efficient_designs`, `efficient_value`, `conditional_*`, `perform_ensemble_designs`) now accept an `rng::AbstractRNG` keyword threaded into `Sim` for reproducibility; the `DistanceBased` sampler tightens its `rng` parameter to `AbstractRNG`. Tests: add Pareto-front edge cases (ties, duplicates, chain, tolerance) and a Mahalanobis insertion-order regression. Update sampler `applicable` checks to pass an RNG instance. Deferred minor findings are recorded in papers/joss/code_review_findings.md. Co-Authored-By: Claude Opus 4.7 --- .../ConditionalUncertaintyReductionMDP.jl | 12 ++- src/GenerativeDesigns/EfficientValueMDP.jl | 9 +- .../UncertaintyReductionMDP.jl | 10 ++- src/GenerativeDesigns/distancebased.jl | 63 ++++++++----- src/StaticDesigns/StaticDesigns.jl | 2 +- src/fronts.jl | 90 ++++++++++++++----- test/GenerativeDesigns/test_distances_sum.jl | 2 +- test/GenerativeDesigns/test_mahalanobis.jl | 13 ++- test/fronts.jl | 15 ++++ 9 files changed, 160 insertions(+), 56 deletions(-) diff --git a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl index 5c11842..9001b8d 100644 --- a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl @@ -132,8 +132,8 @@ end # --- POMDPs interface --- function POMDPs.actions(m::ConditionalUncertaintyReductionMDP, state) all_actions = filter!(collect(keys(m.costs))) do a - !isempty(m.costs[a].features) && - !in(first(m.costs[a].features), keys(state.evidence)) + feats = m.costs[a].features + !isempty(feats) && !any(f -> haskey(state.evidence, f), feats) end return if !isempty(all_actions) && (length(state.evidence) < m.max_experiments) @@ -215,6 +215,7 @@ function conditional_efficient_design( repetitions = 0, realized_uncertainty = false, mdp_options = (;), + rng::AbstractRNG = default_rng(), ) mdp = ConditionalUncertaintyReductionMDP( costs; @@ -245,7 +246,7 @@ function conditional_efficient_design( action, info = action_info(planner, mdp.initial_state) return if repetitions > 0 - queue = [Sim(mdp, planner) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs return (; @@ -305,7 +306,9 @@ function conditional_efficient_designs( repetitions = 0, realized_uncertainty = false, mdp_options = (;), + rng::AbstractRNG = default_rng(), ) + thresholds < 2 && throw(ArgumentError("`thresholds` must be at least 2 (got $thresholds); use `conditional_efficient_design` for a single threshold.")) designs = [] for threshold in range(0.0, 1.0, thresholds) @info "Current threshold level : $threshold" @@ -324,6 +327,7 @@ function conditional_efficient_designs( repetitions, realized_uncertainty, mdp_options, + rng, ), ) end @@ -345,6 +349,7 @@ function perform_ensemble_designs( mdp_options = (;), thred_set = [0.9], N = 30, + rng::AbstractRNG = default_rng(), ) results = Dict() @@ -365,6 +370,7 @@ function perform_ensemble_designs( solver = solver, repetitions = repetitions, mdp_options = mdp_options, + rng = rng, ) push!(ensemble_results, design) end diff --git a/src/GenerativeDesigns/EfficientValueMDP.jl b/src/GenerativeDesigns/EfficientValueMDP.jl index 08135a6..f26d88b 100644 --- a/src/GenerativeDesigns/EfficientValueMDP.jl +++ b/src/GenerativeDesigns/EfficientValueMDP.jl @@ -47,7 +47,7 @@ struct EfficientValueMDP <: POMDPs.MDP{State, Vector{String}} # Check if `sampler`, `uncertainty` are compatible @assert hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) """`sampler` must implement a method accepting `(evidence, readout features, rng)` as its arguments.""" - @assert hasmethod(value, Tuple{Evidence, Vector{Float64}}) """`value` must implement a method accepting `(evidence, costs)` as its argument.""" + @assert hasmethod(value, Tuple{Evidence, NTuple{2, Float64}}) """`value` must implement a method accepting `(evidence, costs)` as its argument, where `costs::NTuple{2, Float64}` is `(monetary cost, execution time)`.""" # actions and their costs costs = Dict{String, ActionCost}( @@ -76,8 +76,8 @@ end function POMDPs.actions(m::EfficientValueMDP, state) all_actions = filter!(collect(keys(m.costs))) do a - return !isempty(m.costs[a].features) && - !in(first(m.costs[a].features), keys(state.evidence)) + feats = m.costs[a].features + return !isempty(feats) && !any(f -> haskey(state.evidence, f), feats) end return collect(powerset(all_actions, 1, m.max_parallel)) @@ -167,6 +167,7 @@ function efficient_value( solver = default_solver, repetitions = 0, mdp_options = (;), + rng::AbstractRNG = default_rng(), ) mdp = EfficientValueMDP(costs; sampler, value, evidence, mdp_options...) @@ -175,7 +176,7 @@ function efficient_value( action, info = action_info(planner, mdp.initial_state) return if repetitions > 0 - queue = [Sim(mdp, planner) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs diff --git a/src/GenerativeDesigns/UncertaintyReductionMDP.jl b/src/GenerativeDesigns/UncertaintyReductionMDP.jl index e79ce6f..f929f91 100644 --- a/src/GenerativeDesigns/UncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/UncertaintyReductionMDP.jl @@ -108,8 +108,8 @@ const eox = "EOX" function POMDPs.actions(m::UncertaintyReductionMDP, state) all_actions = filter!(collect(keys(m.costs))) do a - return !isempty(m.costs[a].features) && - !in(first(m.costs[a].features), keys(state.evidence)) + feats = m.costs[a].features + return !isempty(feats) && !any(f -> haskey(state.evidence, f), feats) end return if !isempty(all_actions) && (length(state.evidence) < m.max_experiments) @@ -215,6 +215,7 @@ function efficient_design( repetitions = 0, realized_uncertainty = false, mdp_options = (;), + rng::AbstractRNG = default_rng(), ) mdp = UncertaintyReductionMDP( costs; @@ -241,7 +242,7 @@ function efficient_design( action, info = action_info(planner, mdp.initial_state) if repetitions > 0 - queue = [Sim(mdp, planner) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs @@ -348,7 +349,9 @@ function efficient_designs( repetitions = 0, realized_uncertainty = false, mdp_options = (;), + rng::AbstractRNG = default_rng(), ) + thresholds < 2 && throw(ArgumentError("`thresholds` must be at least 2 (got $thresholds); use `efficient_design` for a single threshold.")) designs = [] for threshold in range(0.0, 1.0, thresholds) @info "Current threshold level : $threshold" @@ -364,6 +367,7 @@ function efficient_designs( repetitions, realized_uncertainty, mdp_options, + rng, ), ) end diff --git a/src/GenerativeDesigns/distancebased.jl b/src/GenerativeDesigns/distancebased.jl index dc4b0b9..7ce473e 100644 --- a/src/GenerativeDesigns/distancebased.jl +++ b/src/GenerativeDesigns/distancebased.jl @@ -51,6 +51,11 @@ it returns an internal function of `weights` that computes the fraction of varia """ Variance() = function (data; prior) initial = compute_variance(data; weights = prior) + initial > 0 || throw( + ArgumentError( + "`Variance()` requires the target to have non-zero variance under the prior; got $initial. Provide a non-degenerate target column or supply a custom uncertainty functional.", + ), + ) return weights -> (compute_variance(data; weights) / initial) end @@ -69,6 +74,11 @@ function Entropy() return function (labels; prior) @assert elscitype(labels) <: Multiclass "labels must be of `Multiclass` scitype, but `elscitype(labels)=$(elscitype(labels))`" initial = compute_entropy(labels; weights = prior) + initial > 0 || throw( + ArgumentError( + "`Entropy()` requires the target to have non-zero entropy under the prior; got $initial. Provide labels with at least two represented classes or supply a custom uncertainty functional.", + ), + ) return weights -> (compute_entropy(labels; weights) / initial) end end @@ -131,24 +141,23 @@ function SquaredMahalanobisDistance(; diagonal = 0) ) compute_distances = function (evidence::Evidence) - evidence_keys = collect(keys(evidence) ∩ non_targets) + # Canonicalize evidence_keys to the column order used when Λ was + # built; otherwise the dot product silently misaligns axes. + evidence_keys = filter(k -> haskey(evidence, k), non_targets) if length(evidence_keys) == 0 return zeros(nrow(data)) end vec_evidence = map(k -> evidence[k], evidence_keys) + Λ_marginal = Λ[Set(evidence_keys)] factor_p_q = length(non_targets) / length(evidence_keys) distances = map(eachrow(data)) do row # We retrieve the precomputed inverse of the covariance matrix coresponding to the "observed part" # See #12 and, in particular, https://www.jstor.org/stable/3559861 - Λ_marginal = Λ[Set(evidence_keys)] - - factor_p_q * dot( - (vec_evidence - Vector(row[evidence_keys])), - Λ_marginal * (vec_evidence - Vector(row[evidence_keys])), - ) + diff = vec_evidence - Vector(row[evidence_keys]) + factor_p_q * dot(diff, Λ_marginal * diff) end return distances @@ -259,33 +268,43 @@ function DistanceBased( end # Convert distances into probabilistic weights. - compute_weights = function (evidence::Evidence) - similarities = prior .* map(x -> similarity(x), compute_distances(evidence)) - - # Perform hard match on target columns. + apply_masks! = function (sims, evidence) + # Hard match on target columns. for colname in collect(keys(evidence)) ∩ targets - similarities .*= data[!, colname] .== evidence[colname] + sims .*= data[!, colname] .== evidence[colname] end - - # Adjustment based on "column-wise" priors. + # Per-column importance / filter masks supplied by the caller. for colname in keys(evidence) if haskey(weights, colname) - similarities .*= weights[colname] + sims .*= weights[colname] end end + return sims + end - # If all similarities were zero, the `Weights` constructor would error. - if sum(similarities) ≈ 0 - similarities .= 1 - for colname in collect(keys(evidence)) ∩ targets - similarities .*= data[!, colname] .== evidence[colname] - end + compute_weights = function (evidence::Evidence) + similarities = prior .* map(x -> similarity(x), compute_distances(evidence)) + apply_masks!(similarities, evidence) + + # Distance-driven similarities can underflow to zero (e.g. evidence far + # from every historical row, or singular-Σ Mahalanobis subsets). Fall + # back to a uniform prior over rows that still satisfy the user-supplied + # hard constraints; if no rows survive those constraints either, fail + # loudly rather than return a meaningless weight vector. + if iszero(sum(similarities)) + similarities .= 1.0 + apply_masks!(similarities, evidence) + iszero(sum(similarities)) && throw( + ArgumentError( + "No historical rows satisfy the hard constraints implied by the current evidence (target hard-match, `filter_range`, or `importance_weights`). Cannot construct a weight vector.", + ), + ) end return Weights(similarities ./ sum(similarities)) end - sampler = function (evidence::Evidence, columns, rng = default_rng()) + sampler = function (evidence::Evidence, columns, rng::AbstractRNG) observed = data[sample(rng, compute_weights(evidence)), :] return Dict(c => observed[c] for c in columns) diff --git a/src/StaticDesigns/StaticDesigns.jl b/src/StaticDesigns/StaticDesigns.jl index b0f5f5f..4995552 100644 --- a/src/StaticDesigns/StaticDesigns.jl +++ b/src/StaticDesigns/StaticDesigns.jl @@ -102,7 +102,7 @@ function evaluate_experiments( X_ = if !isempty(zero_cost_features) X[!, zero_cost_features] else - DataFrame(; dummy = fill(0.0, nrow(data))) + DataFrame(; dummy = fill(0.0, nrow(X))) end perf_eval = evaluate(model, X_, y; kwargs...) diff --git a/src/fronts.jl b/src/fronts.jl index f1215e2..7638a17 100644 --- a/src/fronts.jl +++ b/src/fronts.jl @@ -2,18 +2,32 @@ front(v) front(f, v; atol1=0, atol2=0) -Construct a Pareto front of `v`. Elements of `v` will be masked by `f` in the computation. +Construct the Pareto front of `v` under joint minimization on the two +objectives `f(x)[1]` and `f(x)[2]`. When `f` is omitted, `identity` is used, +so each element of `v` is taken to be a `(c1, c2)` pair directly. -The first and second (objective) coordinates have to differ by at least `atol1`, `atol2`, respectively, relatively to the latest point on the front. +A point `p` is on the front iff no other point strictly dominates it: there +is no `q ∈ v` with `f(q)[1] ≤ f(p)[1]` and `f(q)[2] ≤ f(p)[2]` and at least +one inequality strict. Points with identical `(c1, c2)` are mutually +non-dominated and are all retained. + +The optional tolerances `atol1`, `atol2` thin the front *after* the +dominance pass: a candidate is kept only if its first coordinate exceeds the +last kept point's by at least `atol1` *and* its second coordinate is at +least `atol2` smaller. With both tolerances `0` (the default), no thinning +occurs and the result is the exact Pareto front. # Examples -```jldocstest -v = [(1,2), (2,3), (2,1)] +```jldoctest +v = [(1, 2), (2, 3), (2, 1)] front(v) # output -[(1, 2), (2, 1)] + +2-element Vector{Tuple{Int64, Int64}}: + (1, 2) + (2, 1) ``` ```jldoctest @@ -22,7 +36,9 @@ front(x -> x[2], v) # output -[(1, (1, 2)), (3, (2, 1))] +2-element Vector{Tuple{Int64, Tuple{Int64, Int64}}}: + (1, (1, 2)) + (3, (2, 1)) ``` ```jldoctest @@ -31,7 +47,9 @@ front(v; atol2 = 0.2) # output -[(1, 2), (3, 1)] +2-element Vector{Tuple{Int64, Float64}}: + (1, 2.0) + (3, 1.0) ``` """ function front end @@ -46,27 +64,57 @@ function front( atol1::Float64 = 0.0, atol2::Float64 = atol1, ) where {F <: Function, T <: AbstractVector} - # dict sort + # Sort lexicographically with a strict-weak comparator. v_sorted = sort( v; - lt = (x, y) -> - (f(x)[1] < f(y)[1] || (f(x)[1] == f(y)[1] && f(x)[2] <= f(y)[2])), + lt = (x, y) -> begin + fx, fy = f(x), f(y) + fx[1] == fy[1] ? fx[2] < fy[2] : fx[1] < fy[1] + end, ) - # check if the second coordinate drops below the second coordinate of the last non-dominated point - ix_front = 1 - ix_current = 2 - while ix_current <= length(v_sorted) - if (f(v_sorted[ix_current])[2] < f(v_sorted[ix_front])[2] - atol2) && - (f(v_sorted[ix_current])[1] > f(v_sorted[ix_front])[1] + atol1) - ix_front = ix_current - ix_current += 1 - else - deleteat!(v_sorted, ix_current) + # Online Pareto front (minimization on both axes). Because `v_sorted` is + # ordered by ascending first coordinate then ascending second, a candidate + # is non-dominated by all previously kept points iff its second coordinate + # is strictly below the running minimum, OR it ties exactly with the most + # recently kept point on both coordinates (mutually non-dominated tie). + pareto = empty(v_sorted) + min_y = Inf + for p in v_sorted + fp = f(p) + if isempty(pareto) + push!(pareto, p) + min_y = fp[2] + elseif fp[2] < min_y + push!(pareto, p) + min_y = fp[2] + elseif fp[2] == min_y + flast = f(pareto[end]) + if flast[1] == fp[1] + push!(pareto, p) + end end end - return v_sorted + # Tolerance thinning: keep a candidate only if it is meaningfully separated + # from the last kept front point on both axes. With atol1 = atol2 = 0 (the + # default), this pass is a no-op. + if atol1 == 0.0 && atol2 == 0.0 + return pareto + end + result = empty(pareto) + for p in pareto + if isempty(result) + push!(result, p) + continue + end + fp = f(p) + flast = f(result[end]) + if (fp[1] - flast[1]) >= atol1 && (flast[2] - fp[2]) >= atol2 + push!(result, p) + end + end + return result end """ diff --git a/test/GenerativeDesigns/test_distances_sum.jl b/test/GenerativeDesigns/test_distances_sum.jl index 013a304..efe8b34 100644 --- a/test/GenerativeDesigns/test_distances_sum.jl +++ b/test/GenerativeDesigns/test_distances_sum.jl @@ -35,7 +35,7 @@ r = DistanceBased( # test signatures using Random: default_rng -@test applicable(sampler, evidence, ["HeartDisease"], default_rng) +@test applicable(sampler, evidence, ["HeartDisease"], default_rng()) @test applicable(uncertainty, evidence) @test applicable(weights, evidence) diff --git a/test/GenerativeDesigns/test_mahalanobis.jl b/test/GenerativeDesigns/test_mahalanobis.jl index 0f66e3c..9f43aa7 100644 --- a/test/GenerativeDesigns/test_mahalanobis.jl +++ b/test/GenerativeDesigns/test_mahalanobis.jl @@ -34,7 +34,7 @@ r = DistanceBased( # test signatures using Random: default_rng -@test applicable(sampler, evidence, ["HeartDisease"], default_rng) +@test applicable(sampler, evidence, ["HeartDisease"], default_rng()) @test applicable(uncertainty, evidence) @test applicable(weights, evidence) @@ -120,3 +120,14 @@ design = efficient_value(experiments; sampler, value, evidence, solver, repetiti design = efficient_value(experiments; sampler, value, evidence, solver); @test design isa Tuple @test !hasproperty(design[2], :stats) + +# Regression: the distance must be invariant to evidence-insertion order. +# Prior to canonicalizing `evidence_keys` against `non_targets`, two +# evidences with the same content but different insertion orders produced +# different Mahalanobis distances when Σ had non-uniform off-diagonals. +let + e1 = Evidence("RestingBP" => 130.0, "Cholesterol" => 250.0, "MaxHR" => 150.0) + e2 = Evidence("MaxHR" => 150.0, "Cholesterol" => 250.0, "RestingBP" => 130.0) + @test weights(e1) ≈ weights(e2) + @test uncertainty(e1) ≈ uncertainty(e2) +end diff --git a/test/fronts.jl b/test/fronts.jl index fc34870..af027b4 100644 --- a/test/fronts.jl +++ b/test/fronts.jl @@ -21,6 +21,21 @@ using CEEDesigns v = [(1, 2), (2, 3), (2, 1)] @test front(v) == [(1, 2), (2, 1)] + + # Ties on the first coordinate: the strictly better-y point wins. + @test front([(1, 5), (1, 1)]) == [(1, 1)] + + # Identical (c1, c2) pairs are mutually non-dominated and both retained. + @test front([(1.0, 2.0), (1.0, 2.0)]) == [(1.0, 2.0), (1.0, 2.0)] + + # A chain of small improvements: every point is on the true front. + v_chain = [(1, 10), (2, 9.5), (3, 9.0), (4, 8.5)] + @test front(v_chain) == v_chain + + # Tolerance thinning keeps only meaningfully separated points: from + # (1,10), (2,9.5) is too close (0.5 < 0.6) and is skipped; (3,9.0) is far + # enough (1.0 ≥ 0.6) and is kept; (4,8.5) is then too close to (3,9.0). + @test front(v_chain; atol1 = 0.0, atol2 = 0.6) == [(1, 10), (3, 9.0)] end @testset "`front(f, v)` tests" begin From dd09bd7b9e67bbedce9792cd2a924370dbf868a1 Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 01:42:25 +0200 Subject: [PATCH 02/11] Merge: StaticDesigns Symbol input + arrangements guard --- src/StaticDesigns/StaticDesigns.jl | 12 +++++-- src/StaticDesigns/arrangements.jl | 27 ++++++++++++++ test/StaticDesigns/test.jl | 58 ++++++++++++++++++++++++++++++ 3 files changed, 95 insertions(+), 2 deletions(-) diff --git a/src/StaticDesigns/StaticDesigns.jl b/src/StaticDesigns/StaticDesigns.jl index 4995552..8e8f9d7 100644 --- a/src/StaticDesigns/StaticDesigns.jl +++ b/src/StaticDesigns/StaticDesigns.jl @@ -42,7 +42,9 @@ Evaluations are run in parallel. # Keyword arguments - `max_cardinality`: maximum cardinality of experimental subsets (defaults to the number of experiments). - - `zero_cost_features`: additional zero-cost features available for each experimental subset (defaults to an empty list). + - `zero_cost_features`: additional zero-cost features available for each experimental subset + (defaults to an empty list). Accepts either a `Vector{String}` (e.g. `["Age", "Sex"]`) or a + `Vector{Symbol}` (e.g. `[:Age, :Sex]`); symbols are normalized to strings internally. - `evaluate_empty_subset`: flag indicating whether to evaluate empty experimental subset. A constant column will be added if `zero_cost_features` is empty (defaults to true). - `return_full_metrics`: flag indicating whether to return full `MLJ.PerformanceEvaluation` metrics. Otherwise return an aggregate "measurement" for the first measure (defaults to false). @@ -70,6 +72,8 @@ function evaluate_experiments( return_full_metrics::Bool = false, kwargs..., ) where {T} + # normalize zero-cost feature names so both Vector{String} and Vector{Symbol} work + zero_cost_features = string.(zero_cost_features) # predictive accuracy scores over subsets of experiments scores = Dict{Set{String}, return_full_metrics ? PerformanceEvaluation : Float64}() # generate all the possible subsets from the set of experiments, with a minimum size of 1 and maximum size of 'max_cardinality' @@ -129,7 +133,9 @@ Evaluations are run in parallel. # Keyword arguments - - `zero_cost_features`: additional zero-cost features available for each experimental subset (defaults to an empty list). + - `zero_cost_features`: additional zero-cost features available for each experimental subset + (defaults to an empty list). Accepts either a `Vector{String}` (e.g. `["Age", "Sex"]`) or a + `Vector{Symbol}` (e.g. `[:Age, :Sex]`); symbols are normalized to strings internally. - `evaluate_empty_subset`: flag indicating whether to evaluate empty experimental subset. # Example @@ -144,6 +150,8 @@ function evaluate_experiments( zero_cost_features = [], evaluate_empty_subset::Bool = true, ) where {T} + # normalize zero-cost feature names so both Vector{String} and Vector{Symbol} work + zero_cost_features = string.(zero_cost_features) scores = Dict{Set{String}, ExperimentalEval}() for exp_set in powerset(collect(keys(experiments)), 1) diff --git a/src/StaticDesigns/arrangements.jl b/src/StaticDesigns/arrangements.jl index 020b52e..3297d7d 100644 --- a/src/StaticDesigns/arrangements.jl +++ b/src/StaticDesigns/arrangements.jl @@ -70,6 +70,33 @@ function optimal_arrangement( tradeoff = (1, 0), mdp_kwargs = default_mdp_kwargs, ) where {T} + # validate that `evals` covers every non-terminal state the MDP rollout may visit. + # the rollout visits all subsets of `experiments` of size 0..|experiments|-1, so each + # such subset must be present in `evals` (its `filtration` is read in `reward`). + # This guards against KeyError when callers pass `evals` produced via `evaluate_experiments` + # with `max_cardinality` smaller than the full experiment set. + missing_states = Set{String}[] + for sub in powerset(collect(experiments), 0, length(experiments) - 1) + s = Set{String}(sub) + if !haskey(evals, s) + push!(missing_states, s) + end + end + if !isempty(missing_states) + n_missing = length(missing_states) + example = first(missing_states) + throw( + ArgumentError( + "`evals` is missing $(n_missing) experimental subset(s) required by the " * + "arrangement search (e.g. $(example)). This typically happens when " * + "`evaluate_experiments` was called with a `max_cardinality` smaller than " * + "the full experiment set. Re-run `evaluate_experiments` without " * + "`max_cardinality` (and with `evaluate_empty_subset = true`) before " * + "computing arrangements over the full set.", + ), + ) + end + experimental_costs = Dict{String, NTuple{2, Float64}}( k => if v isa Number (convert(Float64, v), v) diff --git a/test/StaticDesigns/test.jl b/test/StaticDesigns/test.jl index 2246ba6..ce9f3cd 100644 --- a/test/StaticDesigns/test.jl +++ b/test/StaticDesigns/test.jl @@ -86,3 +86,61 @@ experiments = Dict( "BloodCholesterol" => (1.0, 20.0) => ["Cholesterol"], "BloodSugar" => (1.0, 20.0) => ["FastingBS"], ) + +## Regression test for #11: zero_cost_features as Vector{Symbol} should work. +# Reuse the binary dataset / experiments defined above. +data = CSV.File("StaticDesigns/data/heart_binary.csv") |> DataFrame +experiments_simple = Dict( + "BloodPressure" => 1.0 => ["RestingBP"], + "ECG" => 5.0 => ["RestingECG", "Oldpeak", "ST_Slope", "MaxHR"], + "BloodCholesterol" => 20.0 => ["Cholesterol"], + "BloodSugar" => 20 => ["FastingBS"], +) + +# (1) Symbol vector input — regression test for the MethodError documented in finding #11. +perf_eval_sym = evaluate_experiments( + experiments_simple, + data; + zero_cost_features = [:Age, :Sex, :ChestPainType, :ExerciseAngina], +) +@test perf_eval_sym isa + Dict{Set{String}, NamedTuple{(:loss, :filtration), Tuple{Float64, Float64}}} + +# (2) String vector input — no regression. +perf_eval_str = evaluate_experiments( + experiments_simple, + data; + zero_cost_features = ["Age", "Sex", "ChestPainType", "ExerciseAngina"], +) +@test perf_eval_str isa + Dict{Set{String}, NamedTuple{(:loss, :filtration), Tuple{Float64, Float64}}} +# Both inputs should yield the same evaluation keys/values. +@test keys(perf_eval_sym) == keys(perf_eval_str) +@test all(perf_eval_sym[k] == perf_eval_str[k] for k in keys(perf_eval_sym)) + +## optimal_arrangement golden path with full evals. +seed!(1) +costs_only = Dict{String, Float64}(k => v[1] for (k, v) in experiments_simple) +arr = CEEDesigns.StaticDesigns.optimal_arrangement( + costs_only, + perf_eval_str, + Set(keys(experiments_simple)); + mdp_kwargs = (; n_iterations = 200, depth = 5, exploration_constant = 3.0), +) +@test arr.arrangement isa Vector{Set{String}} +@test !isempty(arr.arrangement) +@test reduce(union, arr.arrangement) == Set(keys(experiments_simple)) +@test arr.monetary_cost > 0 + +## M3 regression: optimal_arrangement with `max_cardinality`-truncated evals +## must raise a clear ArgumentError (not KeyError) when called over the full set. +# Simulate the truncation by dropping evals for subsets of size > 2 from a full eval dict. +truncated_evals = Dict( + k => v for (k, v) in perf_eval_str if length(k) <= 2 +) +@test_throws ArgumentError CEEDesigns.StaticDesigns.optimal_arrangement( + costs_only, + truncated_evals, + Set(keys(experiments_simple)); + mdp_kwargs = (; n_iterations = 50, depth = 5, exploration_constant = 3.0), +) From 5330454efe5a1a9a026c213b75c12896fd43f52f Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 01:43:08 +0200 Subject: [PATCH 03/11] Merge: EfficientValueMDP export + default_solver factory + tests # Conflicts: # src/GenerativeDesigns/EfficientValueMDP.jl --- src/GenerativeDesigns/EfficientValueMDP.jl | 4 +- src/GenerativeDesigns/GenerativeDesigns.jl | 14 +- .../UncertaintyReductionMDP.jl | 25 ++- test/GenerativeDesigns/test.jl | 3 + .../GenerativeDesigns/test_efficient_value.jl | 182 ++++++++++++++++++ 5 files changed, 220 insertions(+), 8 deletions(-) create mode 100644 test/GenerativeDesigns/test_efficient_value.jl diff --git a/src/GenerativeDesigns/EfficientValueMDP.jl b/src/GenerativeDesigns/EfficientValueMDP.jl index f26d88b..072b4e0 100644 --- a/src/GenerativeDesigns/EfficientValueMDP.jl +++ b/src/GenerativeDesigns/EfficientValueMDP.jl @@ -129,7 +129,7 @@ Internally, an instance of the `EfficientValueMDP` structure is created and a su - `sampler`: a function of `(evidence, features, rng)`, in which `evidence` denotes the current experimental evidence, `features` represent the set of features we want to sample from, and `rng` is a random number generator; it returns a dictionary mapping the features to outcomes. - `value`: a function of `(evidence, (monetary costs, execution time))`; it quantifies the utility of experimental evidence. - `evidence=Evidence()`: initial experimental evidence. - - `solver=default_solver`: a POMDPs.jl compatible solver used to solve the decision process. The default solver is [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/). + - `solver=default_solver()`: a POMDPs.jl compatible solver used to solve the decision process. The default solver, returned by the [`default_solver`](@ref) factory, is a fresh [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/) per call. - `repetitions=0`: number of runoffs used to estimate the expected experimental cost. - `mdp_options`: a `NamedTuple` of additional keyword arguments that will be passed to the constructor of [`EfficientValueMDP`](@ref). @@ -164,7 +164,7 @@ function efficient_value( sampler, value, evidence = Evidence(), - solver = default_solver, + solver = default_solver(), repetitions = 0, mdp_options = (;), rng::AbstractRNG = default_rng(), diff --git a/src/GenerativeDesigns/GenerativeDesigns.jl b/src/GenerativeDesigns/GenerativeDesigns.jl index 624fe51..7e4af70 100644 --- a/src/GenerativeDesigns/GenerativeDesigns.jl +++ b/src/GenerativeDesigns/GenerativeDesigns.jl @@ -20,7 +20,7 @@ export conditional_likelihood export Variance, Entropy export Evidence, State export efficient_design, efficient_designs -export efficient_value +export EfficientValueMDP, efficient_value export ConditionalUncertaintyReductionMDP export conditional_efficient_design, conditional_efficient_designs export perform_ensemble_designs @@ -58,7 +58,17 @@ const ActionCost = NamedTuple{(:costs, :features), <:Tuple{NTuple{2, Float64}, V const const_bigM = 1_000_000 -const default_solver = DPWSolver(; n_iterations = 100_000, tree_in_info = true) +""" + default_solver() + +Return a fresh `DPWSolver` with the package's default settings. + +This is a zero-argument factory rather than a shared constant so that each call +site receives an independent solver. MCTS-style solvers carry mutable RNG and +internal state, so sharing a single instance across calls (especially with +parallel runoffs) would compromise reproducibility and concurrency safety. +""" +default_solver() = DPWSolver(; n_iterations = 100_000, tree_in_info = true) # minimize the expected experimental cost while ensuring the uncertainty remains below a specified threshold. include("UncertaintyReductionMDP.jl") diff --git a/src/GenerativeDesigns/UncertaintyReductionMDP.jl b/src/GenerativeDesigns/UncertaintyReductionMDP.jl index f929f91..898067f 100644 --- a/src/GenerativeDesigns/UncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/UncertaintyReductionMDP.jl @@ -7,6 +7,23 @@ In this experimental setup, our objective is to minimize the expected experiment Internally, a state of the decision process is modeled as a tuple `(evidence::Evidence, [total accumulated monetary cost, total accumulated execution time])`. +The per-step reward is the negative marginal cost incurred at that step, i.e. +`-(costs_tradeoff' * Δcosts)` where `Δcosts` is the increment in +`(monetary_cost, execution_time)` between successive states. Reaching a state +where uncertainty is still above `threshold` after `max_experiments` produces a +terminal `-bigM` penalty. Consequently: + + - With `discount = 1.0` (the default), the discounted return telescopes and + the planner maximizes `-(costs_tradeoff' * total_costs)`, i.e. it + minimizes the **total** expected experimental cost (subject to the + threshold and `bigM` penalty). + - With `discount < 1.0`, the planner instead maximizes + `-Σₜ γᵗ * (costs_tradeoff' * Δcostsₜ)`, a γ-weighted sum of marginal + cost increments. Cost incurred *earlier* in the sequence is weighted more + heavily than cost incurred later, so the resulting policy is biased toward + deferring expensive experiments rather than minimizing total cost. Use + `discount = 1.0` if you want the objective to be total expected cost. + # Arguments - `costs`: a dictionary containing pairs `experiment => cost`, where `cost` can either be a scalar cost (modelled as a monetary cost) or a tuple `(monetary cost, execution time)`. @@ -174,7 +191,7 @@ In the uncertainty reduction setup, minimize the expected experimental cost whil - `uncertainty`: a function of `evidence`; it returns the measure of variance or uncertainty about the target variable, conditioned on the experimental evidence acquired so far. - `threshold`: uncertainty threshold. - `evidence=Evidence()`: initial experimental evidence. - - `solver=default_solver`: a POMDPs.jl compatible solver used to solve the decision process. The default solver is [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/). + - `solver=default_solver()`: a POMDPs.jl compatible solver used to solve the decision process. The default solver, returned by the [`default_solver`](@ref) factory, is a fresh [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/) per call. - `repetitions=0`: number of runoffs used to estimate the expected experimental cost. - `mdp_options`: a `NamedTuple` of additional keyword arguments that will be passed to the constructor of [`UncertaintyReductionMDP`](@ref). - `realized_uncertainty=false`: whenever the initial state uncertainty is below the selected threshold, return the actual uncertainty of this state. @@ -211,7 +228,7 @@ function efficient_design( uncertainty, threshold, evidence = Evidence(), - solver = default_solver, + solver = default_solver(), repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -308,7 +325,7 @@ Internally, an instance of the `UncertaintyReductionMDP` structure is created fo - `uncertainty`: a function of `evidence`; it returns the measure of variance or uncertainty about the target variable, conditioned on the experimental evidence acquired so far. - `thresholds`: number of thresholds to consider uniformly in the range between 0 and 1, inclusive. - `evidence=Evidence()`: initial experimental evidence. - - `solver=default_solver`: a POMDPs.jl compatible solver used to solve the decision process. The default solver is [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/). + - `solver=default_solver()`: a POMDPs.jl compatible solver used to solve the decision process. The default solver, returned by the [`default_solver`](@ref) factory, is a fresh [`DPWSolver`](https://juliapomdp.github.io/MCTS.jl/dev/dpw/) per call. - `repetitions=0`: number of runoffs used to estimate the expected experimental cost. - `mdp_options`: a `NamedTuple` of additional keyword arguments that will be passed to the constructor of [`UncertaintyReductionMDP`](@ref). - `realized_uncertainty=false`: whenever the initial state uncertainty is below the selected threshold, return the actual uncertainty of this state. @@ -345,7 +362,7 @@ function efficient_designs( uncertainty, thresholds, evidence = Evidence(), - solver = default_solver, + solver = default_solver(), repetitions = 0, realized_uncertainty = false, mdp_options = (;), diff --git a/test/GenerativeDesigns/test.jl b/test/GenerativeDesigns/test.jl index fe7a85a..9e96023 100644 --- a/test/GenerativeDesigns/test.jl +++ b/test/GenerativeDesigns/test.jl @@ -8,3 +8,6 @@ include("test_mahalanobis.jl") include("test_active_sampling.jl") include("test_conditional_constraints.jl") + +# regression / unit tests for `EfficientValueMDP` +include("test_efficient_value.jl") diff --git a/test/GenerativeDesigns/test_efficient_value.jl b/test/GenerativeDesigns/test_efficient_value.jl new file mode 100644 index 0000000..43ed3e7 --- /dev/null +++ b/test/GenerativeDesigns/test_efficient_value.jl @@ -0,0 +1,182 @@ +# Regression / unit tests for `EfficientValueMDP` and `efficient_value`. +# +# These tests guard against: +# * the previously incorrect `hasmethod(value, Tuple{Evidence, Vector{Float64}})` +# assertion (now `Tuple{Evidence, NTuple{2, Float64}}`); +# * `EfficientValueMDP` not being re-exported from `CEEDesigns.GenerativeDesigns`; +# * regressions in the `repetitions = 0` / `repetitions > 0` and `tree_in_info` +# code paths inside `efficient_value`; +# * RNG reproducibility when the solver is constructed with a seeded RNG. + +using Test +using Random: Xoshiro, AbstractRNG +using CSV, DataFrames +using ScientificTypes +using CEEDesigns, CEEDesigns.GenerativeDesigns + +# ----- Export check (Severe #10) ----- +# `EfficientValueMDP` must be reachable from `using CEEDesigns.GenerativeDesigns`. +@test isdefined(GenerativeDesigns, :EfficientValueMDP) +@test EfficientValueMDP === GenerativeDesigns.EfficientValueMDP + +# ----- Data setup, mirroring `test_mahalanobis.jl` ----- +data = CSV.File("GenerativeDesigns/data/heart_disease.csv") |> DataFrame + +types = Dict( + :MaxHR => Continuous, + :Cholesterol => Continuous, + :ChestPainType => Multiclass, + :Oldpeak => Continuous, + :HeartDisease => Multiclass, + :Age => Continuous, + :ST_Slope => Multiclass, + :RestingECG => Multiclass, + :RestingBP => Continuous, + :Sex => Multiclass, + :FastingBS => Continuous, + :ExerciseAngina => Multiclass, +) +data = coerce(data, types) + +(; sampler, uncertainty, weights) = DistanceBased( + data; + target = "HeartDisease", + uncertainty = Entropy(), + similarity = Exponential(; λ = 5), +) + +experiments = Dict( + "BloodPressure" => 1.0 => ["RestingBP"], + "ECG" => 5.0 => ["RestingECG", "Oldpeak", "ST_Slope", "MaxHR"], + "BloodCholesterol" => 20.0 => ["Cholesterol"], + "BloodSugar" => 20.0 => ["FastingBS"], + "HeartDisease" => 100.0, +) + +evidence = Evidence("Age" => 35, "Sex" => "M") + +# Correct value signature: (Evidence, NTuple{2, Float64}) -> Float64 +value_correct = function (evidence, (monetary_cost, execution_time)) + return (1 - uncertainty(evidence)) - 0.005 * monetary_cost +end + +# Incorrect value signature: (Evidence, Vector{Float64}) -> Float64 +# This was previously *accepted* by the buggy assertion and *rejected* by the +# fixed one. +value_vector = function (evidence::Evidence, costs::Vector{Float64}) + return (1 - uncertainty(evidence)) - 0.005 * costs[1] +end + +# ----- Constructor signature checks (regression for the assertion fix) ----- + +# (1) The correct (Evidence, NTuple{2, Float64}) signature must be accepted. +# The buggy assertion `hasmethod(value, Tuple{Evidence, Vector{Float64}})` +# would have rejected this perfectly valid `value`. +@test ( + EfficientValueMDP(experiments; sampler, value = value_correct, evidence) isa + EfficientValueMDP +) + +# (2) A `value` with the wrong (Evidence, Vector{Float64}) signature must be +# rejected. The buggy assertion would have accepted this even though the +# MDP itself would later pass `state.costs::NTuple{2, Float64}` and crash. +@test_throws AssertionError EfficientValueMDP( + experiments; + sampler, + value = value_vector, + evidence, +) + +# ----- End-to-end happy paths ----- + +# Tiny solver for speed. +small_solver = GenerativeDesigns.DPWSolver(; + n_iterations = 50, + depth = 2, + tree_in_info = true, +) + +# repetitions > 0 with tree_in_info = true +design_with_reps_and_tree = efficient_value( + experiments; + sampler, + value = value_correct, + evidence, + solver = small_solver, + repetitions = 3, +) +@test design_with_reps_and_tree isa Tuple +@test length(design_with_reps_and_tree) == 2 +@test design_with_reps_and_tree[1] isa Real +@test hasproperty(design_with_reps_and_tree[2], :stats) +@test hasproperty(design_with_reps_and_tree[2], :tree) +@test hasproperty(design_with_reps_and_tree[2], :arrangement) + +# repetitions = 0 with tree_in_info = true +design_no_reps_with_tree = efficient_value( + experiments; + sampler, + value = value_correct, + evidence, + solver = small_solver, +) +@test design_no_reps_with_tree isa Tuple +@test length(design_no_reps_with_tree) == 2 +@test design_no_reps_with_tree[1] isa Real +@test !hasproperty(design_no_reps_with_tree[2], :stats) +@test hasproperty(design_no_reps_with_tree[2], :tree) +@test hasproperty(design_no_reps_with_tree[2], :arrangement) + +# repetitions > 0 without tree_in_info +small_solver_no_tree = GenerativeDesigns.DPWSolver(; + n_iterations = 50, + depth = 2, + tree_in_info = false, +) +design_with_reps_no_tree = efficient_value( + experiments; + sampler, + value = value_correct, + evidence, + solver = small_solver_no_tree, + repetitions = 3, +) +@test design_with_reps_no_tree isa Tuple +@test hasproperty(design_with_reps_no_tree[2], :stats) +@test !hasproperty(design_with_reps_no_tree[2], :tree) + +# repetitions = 0 without tree_in_info +design_no_reps_no_tree = efficient_value( + experiments; + sampler, + value = value_correct, + evidence, + solver = small_solver_no_tree, +) +@test design_no_reps_no_tree isa Tuple +@test !hasproperty(design_no_reps_no_tree[2], :stats) +@test !hasproperty(design_no_reps_no_tree[2], :tree) + +# ----- RNG reproducibility ----- +# Two calls that share an identically-seeded `Xoshiro` (one per call, so each +# solver gets its own fresh RNG state) must yield identical objective values. +function run_seeded() + solver = GenerativeDesigns.DPWSolver(; + n_iterations = 50, + depth = 2, + tree_in_info = false, + rng = Xoshiro(42), + ) + return efficient_value( + experiments; + sampler, + value = value_correct, + evidence, + solver, + ) +end + +result_a = run_seeded() +result_b = run_seeded() +@test result_a[1] == result_b[1] +@test result_a[2].arrangement == result_b[2].arrangement From 32724e29bbc6a4abee99148b2e208bb98c42e68d Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 02:09:46 +0200 Subject: [PATCH 04/11] Merge: ConditionalMDP scitype validation + cleanup + tests # Conflicts: # test/GenerativeDesigns/test.jl --- .../ConditionalUncertaintyReductionMDP.jl | 65 ++++- test/GenerativeDesigns/test.jl | 3 + test/GenerativeDesigns/test_conditional.jl | 243 ++++++++++++++++++ tutorials/ConditionalUncertaintyReduction.jl | 2 +- 4 files changed, 303 insertions(+), 10 deletions(-) create mode 100644 test/GenerativeDesigns/test_conditional.jl diff --git a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl index 9001b8d..09e319d 100644 --- a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl @@ -8,6 +8,16 @@ It adds: - data::DataFrame: historical data (rows correspond to weights) - terminal_condition = (target_condition::Dict, tau::Float64): constraint ranges + belief threshold +`target_condition` is a `Dict{String, <:AbstractVector}` (or any pair-style +mapping) whose keys are column names in `data` and whose values are +two-element ranges `[rmin, rmax]`. The membership test is the inclusive +range `rmin <= x <= rmax`, so the keyed columns **must be numeric** (i.e. +`eltype(data[!, colname]) <: Real`). A non-numeric column (e.g. a +`Multiclass` categorical) raises `ArgumentError` at construction. Rows +containing `NaN` in any constraint column are excluded from the conditional +likelihood (because `NaN >= rmin` is `false`); this is documented behavior, +not silent dropping. + Behavior: - The MDP is terminal only when uncertainty <= threshold AND @@ -70,6 +80,30 @@ struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} `uncertainty` must implement a method accepting (evidence). """ + # Validate that every column referenced in `target_condition` is numeric. + # The conditional-likelihood mask uses `>=` / `<=` and is only meaningful + # on columns whose elements support inclusive numeric range comparisons. + target_cond_dict = first(terminal_condition) + for (colname, _) in target_cond_dict + if !(string(colname) in names(data)) + throw( + ArgumentError( + "target_condition references column `$(colname)` which is not present in `data`.", + ), + ) + end + col_eltype = eltype(data[!, string(colname)]) + if !(col_eltype <: Real) + throw( + ArgumentError( + "target_condition column `$(colname)` has eltype `$(col_eltype)`, " * + "but conditional ranges (rmin <= x <= rmax) are only supported on numeric columns. " * + "Coerce the column to a numeric scitype (e.g. Continuous) before constructing the MDP.", + ), + ) + end + end + # Parse costs dict into CEED ActionCost format parsed_costs = Dict{String, ActionCost}( try @@ -184,10 +218,7 @@ function POMDPs.transition(m::ConditionalUncertaintyReductionMDP, state, action_ features = vcat(map(a -> m.costs[a].features, action_set)...) return ImplicitDistribution() do rng - # sampler may return Dict OR (Dict, ...) - obs = m.sampler(state.evidence, features, rng) - observation = obs isa Dict ? obs : first(obs) - + observation = m.sampler(state.evidence, features, rng) return merge(state, observation, (cost_m, cost_t)) end end @@ -211,7 +242,7 @@ function conditional_efficient_design( weights, data, terminal_condition = (Dict(), 0.8), - solver = default_solver, + solver = default_solver(), repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -302,7 +333,7 @@ function conditional_efficient_designs( weights, data, terminal_condition = (Dict(), 0.8), - solver = default_solver, + solver = default_solver(), repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -334,6 +365,22 @@ function conditional_efficient_designs( return front(x -> x[1], designs) end +""" + perform_ensemble_designs(costs; ..., thred_set = [0.9], N = 30) + +Run an ensemble of `N` independent `conditional_efficient_designs` calls per +belief threshold `tau` in `thred_set`, returning a +`Dict{Float64, Vector}` keyed by `tau::Float64`. Each value is the vector of +`N` ensemble outputs (each output is itself the Pareto front returned by +`conditional_efficient_designs`). + +Example access pattern: + +```julia +results = perform_ensemble_designs(experiments; ..., thred_set = [0.6, 0.9]) +runs_at_06 = results[0.6] # Vector of length N +``` +""" function perform_ensemble_designs( costs; sampler, @@ -344,14 +391,14 @@ function perform_ensemble_designs( data, terminal_condition = (Dict(), 0.0), realized_uncertainty = false, - solver = default_solver, + solver = default_solver(), repetitions = 0, mdp_options = (;), thred_set = [0.9], N = 30, rng::AbstractRNG = default_rng(), ) - results = Dict() + results = Dict{Float64, Vector}() for tau in thred_set ensemble_results = [] @@ -374,7 +421,7 @@ function perform_ensemble_designs( ) push!(ensemble_results, design) end - results[:belief => tau] = ensemble_results + results[Float64(tau)] = ensemble_results end return results diff --git a/test/GenerativeDesigns/test.jl b/test/GenerativeDesigns/test.jl index 9e96023..1216537 100644 --- a/test/GenerativeDesigns/test.jl +++ b/test/GenerativeDesigns/test.jl @@ -11,3 +11,6 @@ include("test_conditional_constraints.jl") # regression / unit tests for `EfficientValueMDP` include("test_efficient_value.jl") + +# regression / unit tests for `ConditionalUncertaintyReductionMDP` +include("test_conditional.jl") diff --git a/test/GenerativeDesigns/test_conditional.jl b/test/GenerativeDesigns/test_conditional.jl new file mode 100644 index 0000000..aaa661b --- /dev/null +++ b/test/GenerativeDesigns/test_conditional.jl @@ -0,0 +1,243 @@ +using Test +using CSV, DataFrames +using ScientificTypes +using Random +using Random: Xoshiro +import POMDPs + +using CEEDesigns, CEEDesigns.GenerativeDesigns + +@testset "ConditionalUncertaintyReductionMDP" begin + # ------------------------------------------------------------------ + # Shared fixture data: heart_disease.csv with mixed scitypes. + # ------------------------------------------------------------------ + raw = CSV.File("GenerativeDesigns/data/heart_disease.csv") |> DataFrame + types = Dict( + :MaxHR => Continuous, + :Cholesterol => Continuous, + :Oldpeak => Continuous, + :Age => Continuous, + :RestingBP => Continuous, + :Sex => Multiclass, + :FastingBS => Continuous, + ) + raw_typed = coerce(raw, types) + + # Numeric-only view for the "happy path" tests (mirrors test_mahalanobis.jl). + continuous_cols = + filter(colname -> eltype(raw_typed[!, colname]) == Float64, names(raw_typed)) + data = raw_typed[!, continuous_cols ∪ ["HeartDisease"]] + + (; sampler, uncertainty, weights) = DistanceBased( + data; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = SquaredMahalanobisDistance(; diagonal = 1), + ) + + experiments = Dict( + "BloodPressure" => 1.0 => ["RestingBP"], + "ECG" => 5.0 => ["Oldpeak", "MaxHR"], + "BloodCholesterol" => 20.0 => ["Cholesterol"], + "BloodSugar" => 20.0 => ["FastingBS"], + "HeartDisease" => 100.0, + ) + + # tiny solver to keep the test cheap + tiny_solver = GenerativeDesigns.DPWSolver(; n_iterations = 50, tree_in_info = true) + + # ------------------------------------------------------------------ + # Test 1: severe #13 regression — Multiclass column => ArgumentError + # ------------------------------------------------------------------ + @testset "rejects non-numeric target_condition column (severe #13)" begin + # Build a sampler/weights/uncertainty over the typed dataset which + # still contains the Multiclass `Sex` column. + sub = raw_typed[!, ["Age", "RestingBP", "Sex", "HeartDisease"]] + (; sampler = s2, uncertainty = u2, weights = w2) = DistanceBased( + sub; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + ) + # Should fail at construction with an ArgumentError naming `Sex`. + @test_throws ArgumentError ConditionalUncertaintyReductionMDP( + Dict("BloodPressure" => 1.0 => ["RestingBP"]); + sampler = s2, + uncertainty = u2, + threshold = 0.1, + weights = w2, + data = sub, + terminal_condition = (Dict("Sex" => ("M", "M")), 0.5), + ) + # Re-trigger to inspect the message; ArgumentError should mention "Sex". + err = try + ConditionalUncertaintyReductionMDP( + Dict("BloodPressure" => 1.0 => ["RestingBP"]); + sampler = s2, + uncertainty = u2, + threshold = 0.1, + weights = w2, + data = sub, + terminal_condition = (Dict("Sex" => ("M", "M")), 0.5), + ) + nothing + catch e + e + end + @test err isa ArgumentError + @test occursin("Sex", sprint(showerror, err)) + end + + # ------------------------------------------------------------------ + # Test 2: happy path — conditional_efficient_design returns Tuple shape + # ------------------------------------------------------------------ + @testset "conditional_efficient_design happy path" begin + target_condition = Dict("HeartDisease" => [0.0, 1.0]) + design = conditional_efficient_design( + experiments; + sampler, + uncertainty, + threshold = 0.5, + evidence = Evidence(), + weights, + data, + terminal_condition = (target_condition, 0.0), + solver = tiny_solver, + ) + @test design isa Tuple + @test length(design) == 2 + @test design[1] isa Tuple # (cost, threshold) pair + @test length(design[1]) == 2 + @test hasproperty(design[2], :arrangement) # NamedTuple metadata + end + + # ------------------------------------------------------------------ + # Test 3: thresholds = 1 is rejected (severe #6 regression) + # ------------------------------------------------------------------ + @testset "conditional_efficient_designs rejects thresholds = 1 (severe #6)" begin + target_condition = Dict("HeartDisease" => [0.0, 1.0]) + @test_throws ArgumentError conditional_efficient_designs( + experiments; + sampler, + uncertainty, + thresholds = 1, + evidence = Evidence(), + weights, + data, + terminal_condition = (target_condition, 0.0), + solver = tiny_solver, + ) + end + + # ------------------------------------------------------------------ + # Test 4: perform_ensemble_designs key shape (M7) + # ------------------------------------------------------------------ + @testset "perform_ensemble_designs returns Float64-keyed Dict (M7)" begin + target_condition = Dict("HeartDisease" => [0.0, 1.0]) + thred_set = [0.5, 0.9] + results = perform_ensemble_designs( + experiments; + sampler, + uncertainty, + thresholds = 2, + evidence = Evidence(), + weights, + data, + terminal_condition = (target_condition, 0.0), + solver = tiny_solver, + N = 1, + thred_set = thred_set, + ) + @test results isa Dict + @test Set(keys(results)) == Set(Float64.(thred_set)) + @test all(k -> k isa Float64, keys(results)) + for tau in thred_set + @test haskey(results, Float64(tau)) + @test length(results[Float64(tau)]) == 1 # N = 1 + end + end + + # ------------------------------------------------------------------ + # Test 5: RNG reproducibility — seeded ensemble runs match + # ------------------------------------------------------------------ + @testset "RNG reproducibility with seeded sampler" begin + # Build a *fresh* solver for each run (the default solver carries + # mutable RNG state — see M9), and seed Random.GLOBAL_RNG before + # each call so default RNGs read from the same stream. + target_condition = Dict("HeartDisease" => [0.0, 1.0]) + function _run_once() + Random.seed!(42) + solver = GenerativeDesigns.DPWSolver(; + n_iterations = 50, + tree_in_info = true, + rng = Xoshiro(42), + ) + return perform_ensemble_designs( + experiments; + sampler = sampler, + uncertainty = uncertainty, + thresholds = 2, + evidence = Evidence(), + weights = weights, + data = data, + terminal_condition = (target_condition, 0.0), + solver = solver, + N = 1, + thred_set = [0.5], + ) + end + + results_a = _run_once() + results_b = _run_once() + + # Same keys + @test Set(keys(results_a)) == Set(keys(results_b)) + # Same first-axis costs/thresholds across the two runs + for k in keys(results_a) + front_a = results_a[k][1] + front_b = results_b[k][1] + @test length(front_a) == length(front_b) + for i in eachindex(front_a) + @test front_a[i][1] == front_b[i][1] # the (cost, threshold) Tuple + end + end + end + + # ------------------------------------------------------------------ + # Test 6: empty target_condition reduces to baseline behavior + # ------------------------------------------------------------------ + @testset "empty target_condition: terminal iff uncertainty < threshold" begin + # Force a deterministic uncertainty function so we can drive isterminal. + const_unc = ev -> 0.05 + # constant weights & sampler suffice — the path should not consult them. + const_w = ev -> ones(nrow(data)) ./ nrow(data) + const_sampler = (ev, features, rng) -> Dict(features[1] => 0.0) + + # Empty target_condition + uncertainty (0.05) <= threshold (0.1) => terminal + mdp_terminal = ConditionalUncertaintyReductionMDP( + Dict("BloodPressure" => 1.0 => ["RestingBP"]); + sampler = const_sampler, + uncertainty = const_unc, + threshold = 0.1, + evidence = Evidence(), + weights = const_w, + data = data, + terminal_condition = (Dict(), 0.0), + ) + @test POMDPs.isterminal(mdp_terminal, mdp_terminal.initial_state) == true + + # Same MDP but threshold below uncertainty => not terminal + mdp_nonterminal = ConditionalUncertaintyReductionMDP( + Dict("BloodPressure" => 1.0 => ["RestingBP"]); + sampler = const_sampler, + uncertainty = const_unc, + threshold = 0.01, + evidence = Evidence(), + weights = const_w, + data = data, + terminal_condition = (Dict(), 0.0), + ) + @test POMDPs.isterminal(mdp_nonterminal, mdp_nonterminal.initial_state) == false + end +end diff --git a/tutorials/ConditionalUncertaintyReduction.jl b/tutorials/ConditionalUncertaintyReduction.jl index 9ade548..4704f98 100644 --- a/tutorials/ConditionalUncertaintyReduction.jl +++ b/tutorials/ConditionalUncertaintyReduction.jl @@ -286,7 +286,7 @@ for (idx, evidence) in enumerate(state_init_list) ## Process results for each tau value for tau in taus - runs = ensemble_results[:belief => tau] + runs = ensemble_results[Float64(tau)] df_ensemble = ensemble_to_dataframe(runs) if !haskey(all_ensemble_dfs, tau) From 580bb489e15bf5ed7bb87b079d596edb24f227ae Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 02:10:08 +0200 Subject: [PATCH 05/11] Merge: ensemble fronts label fix + tests + minor guards --- src/StaticDesigns/arrangements.jl | 10 ++ src/ensemble_fronts.jl | 166 ++++++++++-------- .../GenerativeDesigns/test_active_sampling.jl | 4 +- test/runtests.jl | 3 + test/test_ensemble_fronts.jl | 78 ++++++++ 5 files changed, 189 insertions(+), 72 deletions(-) create mode 100644 test/test_ensemble_fronts.jl diff --git a/src/StaticDesigns/arrangements.jl b/src/StaticDesigns/arrangements.jl index 3297d7d..9d94e15 100644 --- a/src/StaticDesigns/arrangements.jl +++ b/src/StaticDesigns/arrangements.jl @@ -116,6 +116,16 @@ function optimal_arrangement( while state != experiments next_action = action(planner, state) + # Guard: the MDP's `actions` always returns non-empty subsets when + # `state != experiments`, so this should be unreachable, but defend + # against an empty action to avoid a confusing `maximum` error and to + # surface the issue if a custom planner ever returns one. + isempty(next_action) && throw( + ArgumentError( + "planner returned an empty action set at state $state; cannot make progress", + ), + ) + push!(arrangement, next_action) monetary_cost += evals[state].filtration * sum(a -> experimental_costs[a][1], next_action) diff --git a/src/ensemble_fronts.jl b/src/ensemble_fronts.jl index 30a9b88..b10cb89 100644 --- a/src/ensemble_fronts.jl +++ b/src/ensemble_fronts.jl @@ -13,7 +13,12 @@ function ensemble_to_dataframe(runs) for design in front threshold = design[1][2] actions = if haskey(design[2], :arrangement) - join(sort(design[2][:arrangement]), ",") + # `arrangement` is a `Vector{Vector{String}}` (each step wraps the + # action set in a vector). Flatten to a `Vector{String}` of action + # names before sorting/joining so the column is a clean + # comma-separated string rather than the literal vector form. + flat_actions = reduce(vcat, design[2][:arrangement]; init = String[]) + join(sort(flat_actions), ",") else "" end @@ -35,53 +40,20 @@ function ensemble_to_dataframe(runs) end """ - plot_ensemble_pareto(df::DataFrame, tau::Float64; show_annotations::Bool=true, normalization::Symbol=:global) - -Create Pareto front visualization for ensemble results with viridis color scheme and horizontal dispersion. - -# Arguments - - - `df::DataFrame`: DataFrame with columns Threshold, Action_Set, Frequency, Average_Utility - - - `tau::Float64`: Belief threshold value for labeling - - `show_annotations::Bool`: Whether to show action set annotations (default: true) - - `normalization::Symbol`: Normalization method for frequencies (default: :global) - - + `:global` - Normalize by max frequency across all points (recommended) - + `:per_threshold` - Normalize within each threshold (can be misleading) - + `:none` - Use raw frequency counts - -# Returns - - - Scatter plot showing Pareto front with: - - + Points with uniform size, colored by probability using viridis scheme - + Horizontal dispersion to show distinct action sets per threshold - + Count annotations showing number of action sets per threshold - + MLASP (Maximum Likelihood Action Sets Path) in red - + Optional annotations for action set changes - -# Notes - - - Uses viridis color scheme for better perceptual uniformity and accessibility - - Horizontal spread indicates number of distinct action sets at each threshold - - Global normalization is recommended as it preserves relative frequencies across thresholds + _compute_dispersion(df_copy) + +Internal helper used by [`plot_ensemble_pareto`](@ref). Given an action-set +DataFrame (already coalesced), compute per-threshold horizontal offsets used to +visually spread points sharing the same threshold and a vector of dispersed +x-coordinates aligned with `eachrow(df_copy)`. + +Returns `(dispersion_offsets, x_coords, action_counts)` where +`dispersion_offsets::Dict{Float64, Vector{Float64}}` is keyed by threshold, +`x_coords::Vector{Float64}` aligns with rows of `df_copy`, and +`action_counts::Dict{Float64, Int}` records how many rows fall under each +threshold. """ -function plot_ensemble_pareto( - df::DataFrame, - tau::Float64; - show_annotations::Bool = true, - normalization::Symbol = :global, - ) - # Calculate normalized frequencies for coloring - df_copy = copy(df) - - # Handle missing action sets - if "Action_Set" in names(df_copy) - df_copy.Action_Set = coalesce.(df_copy.Action_Set, "[\"No Action\"]") - end - - # Calculate horizontal dispersion for each threshold to show distinct action sets +function _compute_dispersion(df_copy::DataFrame) threshold_groups = groupby(df_copy, :Threshold) dispersion_offsets = Dict{Float64, Vector{Float64}}() action_counts = Dict{Float64, Int}() @@ -93,20 +65,14 @@ function plot_ensemble_pareto( n_actions = nrow(group) action_counts[threshold] = n_actions - # Calculate systematic offsets to spread points horizontally if n_actions > 1 - # Create evenly spaced offsets - max_offset = min(0.02 * x_range, x_range / 50) # 2% of x-range or x_range/50, whichever is smaller + max_offset = min(0.02 * x_range, x_range / 50) offsets = collect(range(-max_offset, max_offset; length = n_actions)) - # Sort group by frequency to place higher frequency items in center sorted_indices = sortperm(group.Frequency; rev = true) ordered_offsets = similar(offsets) - # Arrange offsets so highest frequency is in center center_first_order = Int[] - left = 1 - right = n_actions for i in 1:n_actions if i % 2 == 1 push!(center_first_order, div(n_actions + i, 2)) @@ -125,32 +91,27 @@ function plot_ensemble_pareto( end end - # Apply dispersion to x-coordinates - # IMPORTANT: Process in the same order to maintain alignment with norm_frequencies + # Apply dispersion to x-coordinates. The n-th row encountered for a given + # threshold while iterating `eachrow(df_copy)` corresponds to position n in + # `dispersion_offsets[threshold]`, because `groupby` preserves within-group + # row order. An explicit per-threshold counter avoids float `==` lookups + # that would collapse duplicate `Average_Utility` rows onto the same offset + # slot. x_coords = Float64[] threshold_indices = Dict{Float64, Int}() - for (idx, row) in enumerate(eachrow(df_copy)) + for row in eachrow(df_copy) threshold = row.Threshold - # Track which index this is within its threshold group if !haskey(threshold_indices, threshold) threshold_indices[threshold] = 1 else threshold_indices[threshold] += 1 end + point_index = threshold_indices[threshold] - # Get the appropriate offset for this point - group_data = filter(r -> r.Threshold == threshold, df_copy) - point_index = findfirst( - r -> - r.Action_Set == row.Action_Set && - r.Average_Utility == row.Average_Utility && - r.Frequency == row.Frequency, - eachrow(group_data), - ) - - if !isnothing(point_index) && haskey(dispersion_offsets, threshold) + if haskey(dispersion_offsets, threshold) && + point_index <= length(dispersion_offsets[threshold]) offset = dispersion_offsets[threshold][point_index] push!(x_coords, row.Average_Utility + offset) else @@ -158,6 +119,68 @@ function plot_ensemble_pareto( end end + return dispersion_offsets, x_coords, action_counts +end + +""" + plot_ensemble_pareto(df::DataFrame, tau::Float64; show_annotations::Bool=true, normalization::Symbol=:global) + +Create Pareto front visualization for ensemble results with viridis color scheme and horizontal dispersion. + +# Arguments + + - `df::DataFrame`: DataFrame with columns Threshold, Action_Set, Frequency, Average_Utility + + - `tau::Float64`: Belief threshold value for labeling + - `show_annotations::Bool`: Whether to show action set annotations (default: true) + - `normalization::Symbol`: Normalization method for frequencies (default: :global) + + + `:global` - Normalize by max frequency across all points (recommended) + + `:per_threshold` - Normalize within each threshold (can be misleading) + + `:none` - Use raw frequency counts + +# Returns + + - Scatter plot showing Pareto front with: + + + Points with uniform size, colored by probability using viridis scheme + + Horizontal dispersion to show distinct action sets per threshold + + Count annotations showing number of action sets per threshold + + MLASP (Maximum Likelihood Action Sets Path) in red + + Optional annotations for action set changes + +# Notes + + - Uses viridis color scheme for better perceptual uniformity and accessibility + - Horizontal spread indicates number of distinct action sets at each threshold + - Global normalization is recommended as it preserves relative frequencies across thresholds +""" +function plot_ensemble_pareto( + df::DataFrame, + tau::Float64; + show_annotations::Bool = true, + normalization::Symbol = :global, + ) + # Guard against empty input. Several downstream calls (e.g. + # `maximum(values(threshold_counts))`, `maximum(df_copy.Average_Utility)`) + # would otherwise raise an opaque `ArgumentError` from `maximum`/`argmax`. + isempty(df) && throw( + ArgumentError( + "plot_ensemble_pareto requires a non-empty DataFrame (got 0 rows)", + ), + ) + + # Calculate normalized frequencies for coloring + df_copy = copy(df) + + # Handle missing action sets + if "Action_Set" in names(df_copy) + df_copy.Action_Set = coalesce.(df_copy.Action_Set, "[\"No Action\"]") + end + + # Calculate horizontal dispersion for each threshold to show distinct action sets + dispersion_offsets, x_coords, action_counts = _compute_dispersion(df_copy) + # Calculate normalized frequencies based on chosen method norm_frequencies = Float64[] colorbar_label = "Probability" @@ -458,7 +481,10 @@ function plot_ensemble_pareto( # Only annotate if action set changes if current_action_set != prev_action_set - # Clean up the label + # Labels produced by `ensemble_to_dataframe` are already a clean + # comma-separated string. The regex below is retained as a + # defensive no-op for any caller that may still pass labels + # containing bracket/quote characters. label_text = string(current_action_set) label_text = replace(label_text, r"[\[\]'\"\\]" => "") if length(label_text) > 30 diff --git a/test/GenerativeDesigns/test_active_sampling.jl b/test/GenerativeDesigns/test_active_sampling.jl index db56ea6..1022cdb 100644 --- a/test/GenerativeDesigns/test_active_sampling.jl +++ b/test/GenerativeDesigns/test_active_sampling.jl @@ -39,7 +39,7 @@ importance_weights = Dict("B" => x -> 2 <= x <= 7) evidence_1 = Evidence("A" => 5) assigned_weights_1 = weights(evidence_1) - # Check if weights are zero only in the given range + # Check that weights are zero outside the desirable range and positive inside @test all(assigned_weights_1[1:2] .== 0.0) @test all(assigned_weights_1[9:10] .== 0.0) @test all(assigned_weights_1[3:8] .> 0) @@ -48,7 +48,7 @@ importance_weights = Dict("B" => x -> 2 <= x <= 7) assigned_weights_2 = weights(evidence_2) @show assigned_weights_2 - # Check if weights are zero only in the given range + # Check that weights are zero outside the importance-weighted range and positive inside @test all(assigned_weights_2[1] == 0.0) @test all(assigned_weights_2[8:10] .== 0) @test all(assigned_weights_2[2:7] .> 0.0) diff --git a/test/runtests.jl b/test/runtests.jl index dabbc62..ada6905 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,9 @@ using SafeTestsets, BenchmarkTools @time @safetestset "fronts tests" begin include("fronts.jl") end + @time @safetestset "ensemble_fronts tests" begin + include("test_ensemble_fronts.jl") + end @time @safetestset "`StaticDesigns` tests" begin include("StaticDesigns/test.jl") end diff --git a/test/test_ensemble_fronts.jl b/test/test_ensemble_fronts.jl new file mode 100644 index 0000000..799758e --- /dev/null +++ b/test/test_ensemble_fronts.jl @@ -0,0 +1,78 @@ +using Test +using DataFrames +using CEEDesigns +using CEEDesigns: ensemble_to_dataframe, plot_ensemble_pareto, _compute_dispersion + +# Build a synthetic "run" matching the shape produced by the GenerativeDesigns +# MDP wrappers (`efficient_design` / `efficient_value`): +# front ::Vector of designs +# design ::Tuple{Tuple{Utility,Threshold}, NamedTuple{(:arrangement, ...)}} +# arrangement::Vector{Vector{String}} (each step wraps the action set) +function make_design(utility, threshold, arrangement) + return ((utility, threshold), (; arrangement = arrangement)) +end + +@testset "ensemble_to_dataframe label cleanliness (M6)" begin + # One run with a single design whose arrangement matches the documented + # shape `Vector{Vector{String}}`. + arrangement = [["BloodPressure"], ["ECG", "Glucose"]] + front = [make_design(1.5, 0.4, arrangement)] + runs = [front] + + df = ensemble_to_dataframe(runs) + + @test nrow(df) == 1 + @test df.Action_Set[1] == "BloodPressure,ECG,Glucose" + # Defensive: must NOT contain the bracket-laden Vector form. + @test !occursin('[', df.Action_Set[1]) + @test !occursin(']', df.Action_Set[1]) + @test !occursin('"', df.Action_Set[1]) + + # Equivalence: two runs whose actions appear in different order should + # collapse onto the same Action_Set key. + runs2 = [ + [make_design(1.0, 0.4, [["B"], ["A"]])], + [make_design(2.0, 0.4, [["A"], ["B"]])], + ] + df2 = ensemble_to_dataframe(runs2) + @test nrow(df2) == 1 + @test df2.Action_Set[1] == "A,B" + @test df2.Frequency[1] == 2 + @test df2.Average_Utility[1] == 1.5 +end + +@testset "_compute_dispersion distributes duplicate Average_Utility (M11)" begin + # Two rows with identical Average_Utility under the same threshold must NOT + # both be assigned the slot-1 offset. + df = DataFrame(; + Threshold = [0.5, 0.5, 0.5], + Action_Set = ["A", "B", "C"], + Average_Utility = [1.0, 1.0, 2.0], + Frequency = [3, 1, 2], + ) + + dispersion_offsets, x_coords, action_counts = _compute_dispersion(df) + + @test action_counts[0.5] == 3 + @test length(dispersion_offsets[0.5]) == 3 + # All three offset slots are distinct (the helper produces a + # `range(-max_offset, max_offset; length = n_actions)` of length 3). + @test length(unique(dispersion_offsets[0.5])) == 3 + # The x_coords vector must produce three distinct values, even though two + # rows share Average_Utility. This is the regression: the old findfirst + + # float-equality logic collapsed both duplicates onto offset slot 1, so + # x_coords[1] == x_coords[2]. With the per-threshold counter fix they + # consume slots 1 and 2 respectively. + @test length(x_coords) == 3 + @test x_coords[1] != x_coords[2] +end + +@testset "plot_ensemble_pareto guards empty input (M12)" begin + empty_df = DataFrame(; + Threshold = Float64[], + Action_Set = String[], + Average_Utility = Float64[], + Frequency = Int[], + ) + @test_throws ArgumentError plot_ensemble_pareto(empty_df, 0.5) +end From ed2689fd0b0dd2feb13cf07ac234065faf70e496 Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 02:11:28 +0200 Subject: [PATCH 06/11] Merge: DistanceBased narrow catch + Mahalanobis lazy/singular detection # Conflicts: # src/GenerativeDesigns/distancebased.jl # test/GenerativeDesigns/test_mahalanobis.jl --- src/GenerativeDesigns/distancebased.jl | 114 ++++++++++++++++----- test/GenerativeDesigns/test_mahalanobis.jl | 92 +++++++++++++++++ 2 files changed, 179 insertions(+), 27 deletions(-) diff --git a/src/GenerativeDesigns/distancebased.jl b/src/GenerativeDesigns/distancebased.jl index 7ce473e..cc2e1a4 100644 --- a/src/GenerativeDesigns/distancebased.jl +++ b/src/GenerativeDesigns/distancebased.jl @@ -130,19 +130,69 @@ function SquaredMahalanobisDistance(; diagonal = 0) @warn "Not all column types in the predictor matrix are numeric ($(eltype.(eachcol(data)))). This may cause errors." end - Λ = Dict( - Set(features) => begin - Σ = cov(Matrix(data[!, features]), Weights(prior); corrected = false) - # Add diagonal entries. - foreach(i -> Σ[i, i] += diagonal, axes(Σ, 1)) - - inv(Σ) - end for features in powerset(non_targets, 1, length(non_targets)) - ) + # Lazy, memoized cache of inv(Σ_subset). Keys are canonical (sorted-against-`non_targets`) + # vectors of feature names. We avoid the previous eager `2^p − 1` precomputation, and + # we explicitly handle singular submatrices instead of silently returning Inf/NaN + # (Julia's `inv` does not throw on singular matrices). + Λ_cache = Dict{Vector{String}, Matrix{Float64}}() + + get_Λ = function (features::Vector{String}) + # `features` is assumed canonical (sorted against `non_targets`); callers are + # responsible for canonicalization so the cache key round-trips. + cached = get(Λ_cache, features, nothing) + if cached !== nothing + return cached + end + + Σ = cov(Matrix(data[!, features]), Weights(prior); corrected = false) + # Add diagonal entries (regularization, if requested). + foreach(i -> Σ[i, i] += diagonal, axes(Σ, 1)) + + # Detect a singular / near-singular Σ. `inv(Σ)` on a singular matrix may + # silently return Inf/NaN entries in Julia, or throw `SingularException`, + # depending on the factorization path. We handle both by attempting the + # inversion and inspecting the result for non-finite entries. + Λ_marginal = nothing + singular = false + try + Λ_marginal = inv(Σ) + if !all(isfinite, Λ_marginal) + singular = true + end + catch err + if err isa LinearAlgebra.SingularException + singular = true + else + rethrow() + end + end + + if singular + if diagonal > 0 + # User requested regularization but Σ + diagonal·I is still singular: + # this would only happen with extreme inputs; report clearly. + throw( + ArgumentError( + """SquaredMahalanobisDistance: covariance submatrix for features $(features) is singular even after adding `diagonal=$(diagonal)` to the diagonal. Increase `diagonal` or remove collinear/constant columns.""", + ), + ) + else + throw( + ArgumentError( + """SquaredMahalanobisDistance: covariance submatrix for features $(features) is singular (likely due to a collinear or constant column). Pass `diagonal > 0` to `SquaredMahalanobisDistance` for regularization, or remove the offending column from `data`.""", + ), + ) + end + end + + Λ_cache[features] = Λ_marginal + return Λ_marginal + end compute_distances = function (evidence::Evidence) - # Canonicalize evidence_keys to the column order used when Λ was - # built; otherwise the dot product silently misaligns axes. + # Canonicalize evidence_keys against `non_targets` so the cache key, the + # `vec_evidence` ordering, and the row-slice ordering all agree (otherwise + # the dot product silently misaligns axes). evidence_keys = filter(k -> haskey(evidence, k), non_targets) if length(evidence_keys) == 0 @@ -150,12 +200,14 @@ function SquaredMahalanobisDistance(; diagonal = 0) end vec_evidence = map(k -> evidence[k], evidence_keys) - Λ_marginal = Λ[Set(evidence_keys)] + + # Retrieve (and lazily compute) the inverse of the covariance matrix + # corresponding to the "observed part". See #12 and + # https://www.jstor.org/stable/3559861. + Λ_marginal = get_Λ(evidence_keys) factor_p_q = length(non_targets) / length(evidence_keys) distances = map(eachrow(data)) do row - # We retrieve the precomputed inverse of the covariance matrix coresponding to the "observed part" - # See #12 and, in particular, https://www.jstor.org/stable/3559861 diff = vec_evidence - Vector(row[evidence_keys]) factor_p_q * dot(diff, Λ_marginal * diff) end @@ -222,21 +274,29 @@ function DistanceBased( if distance isa Dict distances = Dict( - try - if haskey(distance, colname) - string(colname) => distance[colname] - elseif elscitype(data[!, colname]) <: Continuous - string(colname) => QuadraticDistance() - elseif elscitype(data[!, colname]) <: Multiclass - string(colname) => DiscreteDistance() + begin + if haskey(distance, colname) + # User-supplied distance: do NOT wrap in try/catch — let user code throw + # its own errors so the stacktrace is preserved (see Severe #12). + string(colname) => distance[colname] else - error() + # Only the elscitype-default-dispatch path can legitimately fail with a + # clear "unsupported scitype" message; restrict the try/catch to it. + try + if elscitype(data[!, colname]) <: Continuous + string(colname) => QuadraticDistance() + elseif elscitype(data[!, colname]) <: Multiclass + string(colname) => DiscreteDistance() + else + error() + end + catch + error( + """column $colname has scitype $(elscitype(data[!, colname])), which is not supported by default. + Please provide a custom readout-column distances functional of the signature `(x, col; prior)`.""", + ) + end end - catch - error( - """column $colname has scitype $(elscitype(data[!, colname])), which is not supported by default. - Please provide a custom readout-column distances functional of the signature `(x, col; prior)`.""", - ) end for colname in names(data[!, Not(target)]) ) diff --git a/test/GenerativeDesigns/test_mahalanobis.jl b/test/GenerativeDesigns/test_mahalanobis.jl index 9f43aa7..dbfd5e2 100644 --- a/test/GenerativeDesigns/test_mahalanobis.jl +++ b/test/GenerativeDesigns/test_mahalanobis.jl @@ -131,3 +131,95 @@ let @test weights(e1) ≈ weights(e2) @test uncertainty(e1) ≈ uncertainty(e2) end + +# --- Severe #12: custom `distance` Dict closure must propagate its own errors --- + +struct CustomDistanceError <: Exception end + +@testset "Custom distance closure errors propagate" begin + # User-supplied closure that throws a distinctive error. The previous implementation + # wrapped this branch in a generic `try/catch` that masked the error as "unsupported + # scitype". The error must now surface with its original message. + bad_distance = (x, col; prior = ones(length(col))) -> throw(CustomDistanceError()) + + # The custom-distance branch must NOT be wrapped in a try/catch that re-maps the + # error to the generic "unsupported scitype" message. The original exception type + # has to surface unchanged. + r_bad = DistanceBased( + data; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = Dict("MaxHR" => bad_distance), + ) + @test_throws CustomDistanceError r_bad.weights(Evidence("MaxHR" => 150.0)) +end + +# --- Severe #15: collinear column => singular covariance --- + +@testset "SquaredMahalanobisDistance: collinear column handling" begin + # Synthetic collinear column: exact duplicate of MaxHR. + data_collinear = copy(data) + data_collinear[!, "MaxHR_dup"] = data_collinear[!, "MaxHR"] + + # diagonal == 0 (default) => singular subset Σ; the failure must be explicit. + (; weights) = DistanceBased( + data_collinear; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = SquaredMahalanobisDistance(; diagonal = 0), + ) + # Evidence covering the collinear pair triggers the singular submatrix. + bad_evidence = Evidence("MaxHR" => 150.0, "MaxHR_dup" => 150.0) + @test_throws ArgumentError weights(bad_evidence) + + # diagonal > 0 => regularization kicks in, distances/weights are finite. + (; weights) = DistanceBased( + data_collinear; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = SquaredMahalanobisDistance(; diagonal = 1.0), + ) + ok_weights = weights(bad_evidence) + @test all(isfinite, ok_weights) + @test sum(ok_weights) ≈ 1.0 +end + +# --- Severe #15: lazy cache returns the same numbers as a fresh instance --- + +@testset "SquaredMahalanobisDistance: lazy-cache equivalence" begin + # Build one shared instance whose internal cache will be populated by repeated calls, + # and a second "control" instance for each evidence so we can compare against a freshly + # computed reference (independent cache). + shared = DistanceBased( + data; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = SquaredMahalanobisDistance(; diagonal = 1.0), + ) + + fresh_weights(ev) = DistanceBased( + data; + target = "HeartDisease", + uncertainty = Variance(), + similarity = Exponential(; λ = 5), + distance = SquaredMahalanobisDistance(; diagonal = 1.0), + ).weights(ev) + + ev_a = Evidence("MaxHR" => 150.0) + ev_b = Evidence("MaxHR" => 150.0, "Cholesterol" => 200.0) + + # First pass: populate the cache. + w_a1 = shared.weights(ev_a) + w_b1 = shared.weights(ev_b) + # Second pass on `ev_a`: must hit the cache and return numerically identical weights. + w_a2 = shared.weights(ev_a) + @test w_a1 == w_a2 + + # The cached weights must match a freshly built instance (no shared cache state). + @test w_a1 ≈ fresh_weights(ev_a) + @test w_b1 ≈ fresh_weights(ev_b) +end From df50f5a37dca3d17c6fc8399fc86c19a9ccea11d Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 02:32:22 +0200 Subject: [PATCH 07/11] Fix invalid NamedTuple destructure-and-rename in conditional test Julia's `(; a, b) = nt` syntax does not support `(; a = newname)`. Replace with explicit field accesses for clarity. Co-Authored-By: Claude Opus 4.7 --- test/GenerativeDesigns/test_conditional.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/GenerativeDesigns/test_conditional.jl b/test/GenerativeDesigns/test_conditional.jl index aaa661b..8dd844f 100644 --- a/test/GenerativeDesigns/test_conditional.jl +++ b/test/GenerativeDesigns/test_conditional.jl @@ -54,12 +54,13 @@ using CEEDesigns, CEEDesigns.GenerativeDesigns # Build a sampler/weights/uncertainty over the typed dataset which # still contains the Multiclass `Sex` column. sub = raw_typed[!, ["Age", "RestingBP", "Sex", "HeartDisease"]] - (; sampler = s2, uncertainty = u2, weights = w2) = DistanceBased( + r2 = DistanceBased( sub; target = "HeartDisease", uncertainty = Variance(), similarity = Exponential(; λ = 5), ) + s2, u2, w2 = r2.sampler, r2.uncertainty, r2.weights # Should fail at construction with an ArgumentError naming `Sex`. @test_throws ArgumentError ConditionalUncertaintyReductionMDP( Dict("BloodPressure" => 1.0 => ["RestingBP"]); From 0999e52785cfba2d8bb7610fdce1b1d088a98dfc Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Tue, 26 May 2026 02:38:01 +0200 Subject: [PATCH 08/11] Parameterize MDP structs and drop unused deps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit M2: Replace `::Function` fields on `UncertaintyReductionMDP`, `EfficientValueMDP`, and `ConditionalUncertaintyReductionMDP` with type parameters (e.g. `UncertaintyReductionMDP{S, U}` over the sampler and uncertainty closures). This removes dynamic dispatch on hot MCTS rollout paths. M10: Remove `Reexport` and `Requires` from `Project.toml` deps and `[compat]` — neither is `using`-imported anywhere in `src/`. Co-Authored-By: Claude Opus 4.7 --- Project.toml | 4 ---- .../ConditionalUncertaintyReductionMDP.jl | 18 +++++++++--------- src/GenerativeDesigns/EfficientValueMDP.jl | 14 +++++++------- .../UncertaintyReductionMDP.jl | 14 +++++++------- 4 files changed, 23 insertions(+), 27 deletions(-) diff --git a/Project.toml b/Project.toml index 8ab9363..5a431c4 100644 --- a/Project.toml +++ b/Project.toml @@ -13,8 +13,6 @@ POMDPTools = "7588e00f-9cae-40de-98dc-e0c70c48cdd7" POMDPs = "a93abf59-7444-517b-a68a-c42f96afdd7d" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -Requires = "ae029012-a4dd-5104-9daa-d747884805df" ScientificTypes = "321657f4-b219-11e9-178b-2701a2544e81" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -30,8 +28,6 @@ POMDPTools = "0.1" POMDPs = "0.9" Plots = "1.39" Random = "1.9" -Reexport = "1.2" -Requires = "1.3" ScientificTypes = "3.0" Statistics = "1.9" StatsBase = "0.34" diff --git a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl index 09e319d..bfd08be 100644 --- a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl @@ -25,7 +25,7 @@ Behavior: - Transitions always incorporate sampled evidence; feasibility is enforced through the terminal condition and reward-driven solver behavior. """ -struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} +struct ConditionalUncertaintyReductionMDP{S, U, W} <: POMDPs.MDP{State, Vector{String}} # initial state initial_state::State # uncertainty threshold @@ -45,12 +45,12 @@ struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} bigM::Int64 # sample readouts from the posterior - sampler::Function + sampler::S # measure of uncertainty about the ground truth - uncertainty::Function + uncertainty::U # NEW: compute current state weights - weights::Function + weights::W # NEW: historical data data::DataFrame # NEW: (target constraints, belief threshold tau) @@ -58,8 +58,8 @@ struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} function ConditionalUncertaintyReductionMDP( costs; - sampler, - uncertainty, + sampler::S, + uncertainty::U, threshold, evidence = Evidence(), costs_tradeoff = (1.0, 0.0), @@ -67,10 +67,10 @@ struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} discount = 1.0, bigM = const_bigM, max_experiments = bigM, - weights, + weights::W, data, terminal_condition = (Dict(), 0.0), - ) + ) where {S, U, W} state = State((evidence, Tuple(zeros(2)))) @assert hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) """ @@ -125,7 +125,7 @@ struct ConditionalUncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} end for action in costs ) - return new( + return new{S, U, W}( state, threshold, parsed_costs, diff --git a/src/GenerativeDesigns/EfficientValueMDP.jl b/src/GenerativeDesigns/EfficientValueMDP.jl index 072b4e0..d118b93 100644 --- a/src/GenerativeDesigns/EfficientValueMDP.jl +++ b/src/GenerativeDesigns/EfficientValueMDP.jl @@ -19,7 +19,7 @@ Internally, the reward associated with a particular experimental `evidence` and - `max_parallel`: maximum number of parallel experiments. - `discount`: this is the discounting factor utilized in reward computation. """ -struct EfficientValueMDP <: POMDPs.MDP{State, Vector{String}} +struct EfficientValueMDP{S, V} <: POMDPs.MDP{State, Vector{String}} # initial state initial_state::State @@ -31,18 +31,18 @@ struct EfficientValueMDP <: POMDPs.MDP{State, Vector{String}} discount::Float64 ## sample readouts from the posterior - sampler::Function + sampler::S ## measure of utility - value::Function + value::V function EfficientValueMDP( costs; - sampler, - value, + sampler::S, + value::V, evidence = Evidence(), max_parallel::Int = 1, discount = 1.0, - ) + ) where {S, V} state = State((evidence, Tuple(zeros(2)))) # Check if `sampler`, `uncertainty` are compatible @@ -70,7 +70,7 @@ struct EfficientValueMDP <: POMDPs.MDP{State, Vector{String}} end for action in costs ) - return new(state, costs, max_parallel, discount, sampler, value) + return new{S, V}(state, costs, max_parallel, discount, sampler, value) end end diff --git a/src/GenerativeDesigns/UncertaintyReductionMDP.jl b/src/GenerativeDesigns/UncertaintyReductionMDP.jl index 898067f..077213f 100644 --- a/src/GenerativeDesigns/UncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/UncertaintyReductionMDP.jl @@ -40,7 +40,7 @@ terminal `-bigM` penalty. Consequently: - `bigM`: it refers to the penalty that arises in a scenario where further experimental action is not an option, yet the uncertainty exceeds the allowable limit. - `max_experiments`: this denotes the maximum number of experiments that are permissible to be conducted. """ -struct UncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} +struct UncertaintyReductionMDP{S, U} <: POMDPs.MDP{State, Vector{String}} # initial state initial_state::State # uncertainty threshold @@ -60,14 +60,14 @@ struct UncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} bigM::Int64 ## sample readouts from the posterior - sampler::Function + sampler::S ## measure of uncertainty about the ground truth - uncertainty::Function + uncertainty::U function UncertaintyReductionMDP( costs; - sampler, - uncertainty, + sampler::S, + uncertainty::U, threshold, evidence = Evidence(), costs_tradeoff = (1, 0), @@ -75,7 +75,7 @@ struct UncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} discount = 1.0, bigM = const_bigM, max_experiments = bigM, - ) + ) where {S, U} state = State((evidence, Tuple(zeros(2)))) # check if `sampler`, `uncertainty` are compatible @@ -103,7 +103,7 @@ struct UncertaintyReductionMDP <: POMDPs.MDP{State, Vector{String}} end for action in costs ) - return new( + return new{S, U}( state, threshold, costs, From 58519791df55583c1cf8649144871255a35f9ba0 Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Wed, 27 May 2026 17:05:53 +0200 Subject: [PATCH 09/11] Format with Runic --- src/GenerativeDesigns/distancebased.jl | 32 +++++++++++++------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/GenerativeDesigns/distancebased.jl b/src/GenerativeDesigns/distancebased.jl index cc2e1a4..5f4077e 100644 --- a/src/GenerativeDesigns/distancebased.jl +++ b/src/GenerativeDesigns/distancebased.jl @@ -275,29 +275,29 @@ function DistanceBased( if distance isa Dict distances = Dict( begin - if haskey(distance, colname) - # User-supplied distance: do NOT wrap in try/catch — let user code throw - # its own errors so the stacktrace is preserved (see Severe #12). - string(colname) => distance[colname] + if haskey(distance, colname) + # User-supplied distance: do NOT wrap in try/catch — let user code throw + # its own errors so the stacktrace is preserved (see Severe #12). + string(colname) => distance[colname] else - # Only the elscitype-default-dispatch path can legitimately fail with a - # clear "unsupported scitype" message; restrict the try/catch to it. - try - if elscitype(data[!, colname]) <: Continuous - string(colname) => QuadraticDistance() + # Only the elscitype-default-dispatch path can legitimately fail with a + # clear "unsupported scitype" message; restrict the try/catch to it. + try + if elscitype(data[!, colname]) <: Continuous + string(colname) => QuadraticDistance() elseif elscitype(data[!, colname]) <: Multiclass - string(colname) => DiscreteDistance() + string(colname) => DiscreteDistance() else - error() + error() end catch - error( - """column $colname has scitype $(elscitype(data[!, colname])), which is not supported by default. - Please provide a custom readout-column distances functional of the signature `(x, col; prior)`.""", - ) + error( + """column $colname has scitype $(elscitype(data[!, colname])), which is not supported by default. + Please provide a custom readout-column distances functional of the signature `(x, col; prior)`.""", + ) end end - end for colname in names(data[!, Not(target)]) + end for colname in names(data[!, Not(target)]) ) compute_distances = sum_of_distances(data, targets, distances; prior) From 494cdb7068ef26b6b0c9b9eee92a3a8569bb9d26 Mon Sep 17 00:00:00 2001 From: Jan Bima Date: Wed, 27 May 2026 17:42:36 +0200 Subject: [PATCH 10/11] Document default_solver in API reference Resolves the @ref to default_solver added by the EfficientValueMDP / default-solver factory change; without this, Documenter's cross-reference pass errors out (`Cannot resolve @ref for \[`default_solver`\](@ref)`). --- docs/src/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/api.md b/docs/src/api.md index 05469e2..93522a0 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -25,6 +25,7 @@ CEEDesigns.GenerativeDesigns.Entropy CEEDesigns.GenerativeDesigns.efficient_design CEEDesigns.GenerativeDesigns.efficient_designs CEEDesigns.GenerativeDesigns.efficient_value +CEEDesigns.GenerativeDesigns.default_solver ``` ### Distance-Based Sampling From 9daf77fa5401a5f21b9641b1be295f9bd99cd226 Mon Sep 17 00:00:00 2001 From: "Chen, Tianchi" Date: Wed, 3 Jun 2026 23:43:04 -0400 Subject: [PATCH 11/11] Fix correctness, reproducibility, and API issues from pre-JOSS review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Correctness & crashes: - S1: DiscreteDistance was inverted (match -> large distance). Now `y == x ? 0.0 : λ`, consistent with QuadraticDistance, so rows matching the evidence receive the highest similarity weight. - C1: ConditionalUncertaintyReductionMDP normalizes target_condition keys to String (Symbol keys constructed fine but crashed at runtime) and validates ranges; conditional_likelihood is now robust to String/Symbol keys. - C2: accept Union{Missing,<:Real} constraint columns (nonmissingtype check) and exclude missing rows via coalesce(..., false). - M1: max_experiments now counts completed experiments, not raw evidence entries, so prior evidence / multi-feature experiments are not mis-counted. - M4: QuadraticDistance raises a clear ArgumentError on a zero-variance (constant) column and recomputes σ per call (drops cross-column memoization). - N3: DistanceBased validates importance_weights vector length. Reproducibility & concurrency: - S2: default_solver(rng) is seedable; each Sim runoff gets an independent rng (fixes a parallel data race); efficient_designs / conditional_efficient_designs use a fresh solver+rng per threshold; perform_ensemble_designs uses an independent rng per ensemble run. - S3: optimal_arrangement takes an rng (no implicit GLOBAL_RNG); static efficient_designs pre-seeds one rng per task before @threads. Docs / robustness / naming: - C3: public-input @assert -> ArgumentError across MDP constructors, conditional_likelihood, and Entropy. - D1: EfficientValueMDP docstring value signature -> value(evidence, costs). - D2: Exponential docstring (λ=1/2, exp(-λx)); typos; unified terminal_condition default; tau_set alias (thred_set still accepted). - D3: static evaluate_experiments normalizes Symbol features; plot_front empty guard and get_labels -> make_labels docstring fix; ensemble column Average_Utility -> Average_Cost. - M2: ensemble_to_dataframe de-duplicates identical points within a run so realized_uncertainty collapses no longer inflate MLASP frequencies. - P1/P2: documented normalized uncertainty and hard terminal condition. Also adds CHANGELOG.md summarizing these changes, regression tests (test_review_fixes.jl, ensemble M2 dedupe, static S3 reproducibility), and .gitignore entries for local workspace/tooling artifacts. Full suite green with JULIA_NUM_THREADS=4. --- .gitignore | 18 ++- CHANGELOG.md | 99 +++++++++++++++ .../ConditionalUncertaintyReductionMDP.jl | 115 +++++++++++------ src/GenerativeDesigns/EfficientValueMDP.jl | 26 ++-- src/GenerativeDesigns/GenerativeDesigns.jl | 10 +- .../UncertaintyReductionMDP.jl | 48 ++++++-- src/GenerativeDesigns/distancebased.jl | 49 +++++--- src/StaticDesigns/StaticDesigns.jl | 15 ++- src/StaticDesigns/arrangements.jl | 3 +- src/ensemble_fronts.jl | 31 +++-- src/fronts.jl | 5 +- test/GenerativeDesigns/test.jl | 3 + .../GenerativeDesigns/test_efficient_value.jl | 2 +- test/GenerativeDesigns/test_review_fixes.jl | 116 ++++++++++++++++++ test/StaticDesigns/test.jl | 26 ++++ test/test_ensemble_fronts.jl | 28 ++++- 16 files changed, 497 insertions(+), 97 deletions(-) create mode 100644 CHANGELOG.md create mode 100644 test/GenerativeDesigns/test_review_fixes.jl diff --git a/.gitignore b/.gitignore index b9d59cc..1586f00 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,20 @@ LocalPreferences.toml _* tutorials/build .env -.ipynb_checkpoints \ No newline at end of file +.ipynb_checkpoints + +# Local workspace and tooling artifacts (do not commit) +.claude-plugin/ +.craft-agent/ +.agents/ +.pi-agent/ +.pi-sessions/ +/config.json +/events.jsonl +/views.json +/theme.json +/labels/ +/statuses/ +/sessions/ +/skills/ +/sources/ \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..2f36fb5 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,99 @@ +# Changelog + +All notable changes to CEEDesigns.jl are documented in this file. + +## [Unreleased] + +This update resolves a set of correctness, reproducibility, and API/usability +issues identified during a pre-JOSS code review. Each item lists the previous +behavior and the change. Issue tags (S#, C#, M#, N#, D#, P#) are internal +review references. + +### Fixed — correctness & crashes + +- **Categorical distance direction (S1).** `DiscreteDistance` returned `λ` for a + *match* and `0` for a *mismatch* — the opposite of a distance. After collation + through the `Exponential` similarity, historical rows that *matched* the + evidence received the *lowest* weight. It now returns `0` on a match and `λ` on + a mismatch, consistent with `QuadraticDistance`, so matching rows receive the + highest similarity weight. + *Note:* this changes the numerical output of any model that uses categorical + features. + +- **Symbol-keyed constraints crashed at run time (C1).** + `ConditionalUncertaintyReductionMDP` validated `target_condition` keys by + stringifying them but stored the original keys, so a `Dict(:Col => …)` + constructed successfully and then threw "Column … not found" during + `isterminal`. Keys are now normalized to `String` (and ranges validated) at + construction, and `conditional_likelihood` accepts both `String` and `Symbol` + keys. + +- **Columns with missing values were rejected (C2).** The numeric check + `eltype(col) <: Real` excluded `Union{Missing, Float64}` columns. It now uses + `nonmissingtype(...) <: Real`, and the membership mask uses + `coalesce(…, false)` so rows with `missing` in a constraint column are excluded + rather than erroring. + +- **`max_experiments` counted evidence entries (M1).** The budget check counted + evidence dictionary entries (including prior evidence and each feature of a + multi-feature experiment) instead of completed experiments. It now counts + experiments whose features are all present. + +- **Constant feature column produced an opaque error (M4).** A zero-variance + column made `QuadraticDistance` divide by zero, surfacing as the unrelated + "weights cannot contain Inf or NaN values". It now raises a clear + `ArgumentError`, and `σ` is recomputed per call (removing a latent + cross-column memoization bug). + +- **`importance_weights` length was unchecked (N3).** A wrong-length vector failed + later with an opaque broadcasting error. The length is now validated against + `nrow(data)` with a clear `ArgumentError`. + +### Fixed — reproducibility & concurrency + +- **Shared RNG across parallel simulations (S2).** All `Sim` runoffs shared a + single RNG object (a data race under `run_parallel`), and one solver instance + was reused across thresholds and ensemble runs. Each `Sim` now gets an + independent RNG derived from the caller's RNG; `default_solver(rng)` is + seedable; `efficient_designs` / `conditional_efficient_designs` use a fresh + solver and RNG per threshold; and `perform_ensemble_designs` uses an + independent RNG per ensemble run. Results are reproducible from a single seed. + +- **Static arrangement search used the global RNG (S3).** `optimal_arrangement` + built its `DPWSolver` with the implicit `Random.GLOBAL_RNG`, which the threaded + `StaticDesigns.efficient_designs` accessed concurrently. `optimal_arrangement` + now takes an `rng`, and the threaded driver pre-seeds one independent RNG per + task. + +### Changed — API, documentation, robustness + +- **Input validation uses `ArgumentError` (C3).** User-facing `@assert`s in the + MDP constructors, `conditional_likelihood`, and `Entropy` now throw + `ArgumentError` with stable messages. +- **`EfficientValueMDP` value signature (D1).** Docstring corrected from + `value(evidence)` to `value(evidence, costs)`. +- **Documentation / defaults / naming (D2).** `Exponential` docstring corrected to + `λ = 1/2` and `exp(-λ x)`; fixed typos; unified the `terminal_condition` + default; added `tau_set` as the keyword for `perform_ensemble_designs` (the + previous `thred_set` is still accepted). +- **Static designs & plotting (D3).** Experiment feature lists are normalized to + `String` so `Symbol` features work; `plot_front` guards against empty input; a + docstring referenced a nonexistent `get_labels` (now `make_labels`); and the + ensemble output column `Average_Utility` was renamed to `Average_Cost` to + reflect that it holds cost. +- **Ensemble frequency counting (M2).** `ensemble_to_dataframe` now de-duplicates + identical points within a single run, so collapsed points (under + `realized_uncertainty = true`) no longer inflate frequencies / MLASP votes. +- **Documentation note (P1/P2).** Documented that the uncertainty measure is + normalized (thresholds on `[0, 1]`) and that the terminal information term is + realized as a hard terminal condition rather than an additive penalty. + +### Added + +- Regression tests covering the above: `test/GenerativeDesigns/test_review_fixes.jl` + (S1/C1/C2/M1/M4/N3), ensemble within-run de-duplication, and static arrangement + reproducibility. + +### Testing + +- Full test suite passes with `JULIA_NUM_THREADS=4`. diff --git a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl index bfd08be..7d6d4c2 100644 --- a/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/ConditionalUncertaintyReductionMDP.jl @@ -73,36 +73,57 @@ struct ConditionalUncertaintyReductionMDP{S, U, W} <: POMDPs.MDP{State, Vector{S ) where {S, U, W} state = State((evidence, Tuple(zeros(2)))) - @assert hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) """ - `sampler` must implement a method accepting (evidence, features, rng). - """ - @assert hasmethod(uncertainty, Tuple{Evidence}) """ - `uncertainty` must implement a method accepting (evidence). - """ - - # Validate that every column referenced in `target_condition` is numeric. - # The conditional-likelihood mask uses `>=` / `<=` and is only meaningful - # on columns whose elements support inclusive numeric range comparisons. - target_cond_dict = first(terminal_condition) - for (colname, _) in target_cond_dict - if !(string(colname) in names(data)) + hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) || throw( + ArgumentError("`sampler` must implement a method accepting (evidence, features, rng)."), + ) + hasmethod(uncertainty, Tuple{Evidence}) || throw( + ArgumentError("`uncertainty` must implement a method accepting (evidence)."), + ) + + # Validate and normalize `target_condition`. Keys are stringified so Symbol + # column names round-trip to the String column names in `data` (otherwise the + # likelihood lookup fails at runtime). Columns must be numeric (allowing + # `Union{Missing,<:Real}`) because the membership mask uses inclusive + # `rmin <= x <= rmax`. + raw_target_cond = first(terminal_condition) + tau_value = last(terminal_condition) + normalized_target_cond = Dict{String, Tuple{Float64, Float64}}() + for (colname, range) in raw_target_cond + scol = string(colname) + if !(scol in names(data)) throw( ArgumentError( "target_condition references column `$(colname)` which is not present in `data`.", ), ) end - col_eltype = eltype(data[!, string(colname)]) - if !(col_eltype <: Real) + col_eltype = eltype(data[!, scol]) + if !(nonmissingtype(col_eltype) <: Real) throw( ArgumentError( "target_condition column `$(colname)` has eltype `$(col_eltype)`, " * - "but conditional ranges (rmin <= x <= rmax) are only supported on numeric columns. " * + "but conditional ranges (rmin <= x <= rmax) are only supported on numeric columns " * + "(`Union{Missing,<:Real}` is allowed; missing rows are excluded). " * "Coerce the column to a numeric scitype (e.g. Continuous) before constructing the MDP.", ), ) end + if length(range) != 2 + throw( + ArgumentError( + "target_condition for column `$(colname)` must be a 2-element `[rmin, rmax]`; got `$(range)`.", + ), + ) + end + rmin, rmax = Float64(range[1]), Float64(range[2]) + rmin <= rmax || throw( + ArgumentError( + "target_condition for column `$(colname)` has rmin=$(rmin) > rmax=$(rmax).", + ), + ) + normalized_target_cond[scol] = (rmin, rmax) end + terminal_condition = (normalized_target_cond, Float64(tau_value)) # Parse costs dict into CEED ActionCost format parsed_costs = Dict{String, ActionCost}( @@ -151,13 +172,21 @@ function conditional_likelihood( target_condition::Dict, ) w = compute_weights(evidence) - @assert length(w) == nrow(hist_data) "weights length must match number of rows in hist_data" + length(w) == nrow(hist_data) || throw( + ArgumentError( + "weights length ($(length(w))) must match number of rows in hist_data ($(nrow(hist_data))).", + ), + ) valid = trues(nrow(hist_data)) for (colname, range) in target_condition - @assert colname in names(hist_data) "Column $colname not found in data." + scol = string(colname) + scol in names(hist_data) || + throw(ArgumentError("target_condition column `$(colname)` not found in data.")) rmin, rmax = range - valid .&= (hist_data[!, colname] .>= rmin) .& (hist_data[!, colname] .<= rmax) + col = hist_data[!, scol] + # `coalesce(..., false)` excludes rows whose constraint column is `missing`. + valid .&= coalesce.((col .>= rmin) .& (col .<= rmax), false) end return sum(w[valid]) @@ -170,7 +199,15 @@ function POMDPs.actions(m::ConditionalUncertaintyReductionMDP, state) !isempty(feats) && !any(f -> haskey(state.evidence, f), feats) end - return if !isempty(all_actions) && (length(state.evidence) < m.max_experiments) + # Count *completed experiments* (all of an experiment's features present), not + # raw evidence entries, so prior evidence / multi-feature experiments don't + # mis-count against `max_experiments`. + n_completed = count(keys(m.costs)) do a + feats = m.costs[a].features + !isempty(feats) && all(f -> haskey(state.evidence, f), feats) + end + + return if !isempty(all_actions) && (n_completed < m.max_experiments) collect(powerset(all_actions, 1, m.max_parallel)) else [[eox]] @@ -241,8 +278,8 @@ function conditional_efficient_design( evidence = Evidence(), weights, data, - terminal_condition = (Dict(), 0.8), - solver = default_solver(), + terminal_condition = (Dict(), 0.0), + solver = nothing, repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -273,11 +310,11 @@ function conditional_efficient_design( ) end - planner = solve(solver, mdp) + planner = solve(isnothing(solver) ? default_solver(rng) : solver, mdp) action, info = action_info(planner, mdp.initial_state) return if repetitions > 0 - queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng = Xoshiro(rand(rng, UInt64))) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs return (; @@ -332,8 +369,8 @@ function conditional_efficient_designs( evidence = Evidence(), weights, data, - terminal_condition = (Dict(), 0.8), - solver = default_solver(), + terminal_condition = (Dict(), 0.0), + solver = nothing, repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -343,6 +380,8 @@ function conditional_efficient_designs( designs = [] for threshold in range(0.0, 1.0, thresholds) @info "Current threshold level : $threshold" + inner_rng = Xoshiro(rand(rng, UInt64)) + inner_solver = isnothing(solver) ? default_solver(inner_rng) : solver push!( designs, conditional_efficient_design( @@ -354,11 +393,11 @@ function conditional_efficient_designs( weights, data, terminal_condition = terminal_condition, - solver, + solver = inner_solver, repetitions, realized_uncertainty, mdp_options, - rng, + rng = inner_rng, ), ) end @@ -366,10 +405,11 @@ function conditional_efficient_designs( end """ - perform_ensemble_designs(costs; ..., thred_set = [0.9], N = 30) + perform_ensemble_designs(costs; ..., tau_set = [0.9], N = 30) Run an ensemble of `N` independent `conditional_efficient_designs` calls per -belief threshold `tau` in `thred_set`, returning a +belief threshold `tau` in `tau_set` (the legacy name `thred_set` is still +accepted as an alias), returning a `Dict{Float64, Vector}` keyed by `tau::Float64`. Each value is the vector of `N` ensemble outputs (each output is itself the Pareto front returned by `conditional_efficient_designs`). @@ -377,7 +417,7 @@ belief threshold `tau` in `thred_set`, returning a Example access pattern: ```julia -results = perform_ensemble_designs(experiments; ..., thred_set = [0.6, 0.9]) +results = perform_ensemble_designs(experiments; ..., tau_set = [0.6, 0.9]) runs_at_06 = results[0.6] # Vector of length N ``` """ @@ -391,19 +431,24 @@ function perform_ensemble_designs( data, terminal_condition = (Dict(), 0.0), realized_uncertainty = false, - solver = default_solver(), + solver = nothing, repetitions = 0, mdp_options = (;), - thred_set = [0.9], + tau_set = [0.9], + thred_set = nothing, N = 30, rng::AbstractRNG = default_rng(), ) + # `thred_set` is the legacy spelling; prefer `tau_set`. + tau_values = isnothing(thred_set) ? tau_set : thred_set results = Dict{Float64, Vector}() - for tau in thred_set + for tau in tau_values ensemble_results = [] for i in 1:N @info "Running ensemble $i for belief threshold τ=$tau" + # Independent, reproducible rng per ensemble run (derived from `rng`). + run_rng = Xoshiro(rand(rng, UInt64)) design = conditional_efficient_designs( costs; sampler = sampler, @@ -417,7 +462,7 @@ function perform_ensemble_designs( solver = solver, repetitions = repetitions, mdp_options = mdp_options, - rng = rng, + rng = run_rng, ) push!(ensemble_results, design) end diff --git a/src/GenerativeDesigns/EfficientValueMDP.jl b/src/GenerativeDesigns/EfficientValueMDP.jl index d118b93..5efba67 100644 --- a/src/GenerativeDesigns/EfficientValueMDP.jl +++ b/src/GenerativeDesigns/EfficientValueMDP.jl @@ -5,7 +5,7 @@ Structure that parametrizes the experimental decision-making process. It is used In this experimental setup, our objective is to maximize the value of the experimental evidence (such as clinical utility), adjusted for experimental costs. -Internally, the reward associated with a particular experimental `evidence` and with total accumulated `monetary_cost` and (optionally) `execution_time` is computed as `value(evidence) - costs_tradeoff' * [monetary_cost, execution_time]`. +Internally, the reward associated with a particular experimental `evidence` and with total accumulated `monetary_cost` and (optionally) `execution_time` is computed as `value(evidence, costs) - costs_tradeoff' * [monetary_cost, execution_time]`. # Arguments @@ -14,7 +14,7 @@ Internally, the reward associated with a particular experimental `evidence` and # Keyword Arguments - `sampler`: a function of `(evidence, features, rng)`, in which `evidence` denotes the current experimental evidence, `features` represent the set of features we want to sample from, and `rng` is a random number generator; it returns a dictionary mapping the features to outcomes. - - `value`: a function of `(evidence)`; it quantifies the utility of experimental evidence. + - `value`: a function of `(evidence, costs)`, where `costs::NTuple{2, Float64}` is `(monetary cost, execution time)`; it quantifies the utility of experimental evidence. - `evidence=Evidence()`: initial experimental evidence. - `max_parallel`: maximum number of parallel experiments. - `discount`: this is the discounting factor utilized in reward computation. @@ -45,9 +45,17 @@ struct EfficientValueMDP{S, V} <: POMDPs.MDP{State, Vector{String}} ) where {S, V} state = State((evidence, Tuple(zeros(2)))) - # Check if `sampler`, `uncertainty` are compatible - @assert hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) """`sampler` must implement a method accepting `(evidence, readout features, rng)` as its arguments.""" - @assert hasmethod(value, Tuple{Evidence, NTuple{2, Float64}}) """`value` must implement a method accepting `(evidence, costs)` as its argument, where `costs::NTuple{2, Float64}` is `(monetary cost, execution time)`.""" + # Check if `sampler`, `value` are compatible + hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) || throw( + ArgumentError( + "`sampler` must implement a method accepting `(evidence, readout features, rng)` as its arguments.", + ), + ) + hasmethod(value, Tuple{Evidence, NTuple{2, Float64}}) || throw( + ArgumentError( + "`value` must implement a method accepting `(evidence, costs)` as its argument, where `costs::NTuple{2, Float64}` is `(monetary cost, execution time)`.", + ), + ) # actions and their costs costs = Dict{String, ActionCost}( @@ -164,19 +172,19 @@ function efficient_value( sampler, value, evidence = Evidence(), - solver = default_solver(), + solver = nothing, repetitions = 0, mdp_options = (;), rng::AbstractRNG = default_rng(), ) mdp = EfficientValueMDP(costs; sampler, value, evidence, mdp_options...) - # planner - planner = solve(solver, mdp) + # planner (fresh solver seeded by `rng` unless the caller supplied one) + planner = solve(isnothing(solver) ? default_solver(rng) : solver, mdp) action, info = action_info(planner, mdp.initial_state) return if repetitions > 0 - queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng = Xoshiro(rand(rng, UInt64))) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs diff --git a/src/GenerativeDesigns/GenerativeDesigns.jl b/src/GenerativeDesigns/GenerativeDesigns.jl index 7e4af70..7ee4248 100644 --- a/src/GenerativeDesigns/GenerativeDesigns.jl +++ b/src/GenerativeDesigns/GenerativeDesigns.jl @@ -8,7 +8,7 @@ using DataFrames, ScientificTypes using LinearAlgebra using Statistics using StatsBase: Weights, countmap, entropy, sample -using Random: default_rng, AbstractRNG +using Random: default_rng, AbstractRNG, Xoshiro using MCTS using ..CEEDesigns: front @@ -61,14 +61,16 @@ const const_bigM = 1_000_000 """ default_solver() -Return a fresh `DPWSolver` with the package's default settings. +Return a fresh `DPWSolver` with the package's default settings, seeded by `rng`. -This is a zero-argument factory rather than a shared constant so that each call +This is a zero/one-argument factory rather than a shared constant so that each call site receives an independent solver. MCTS-style solvers carry mutable RNG and internal state, so sharing a single instance across calls (especially with parallel runoffs) would compromise reproducibility and concurrency safety. +Passing an `rng` makes the planning step reproducible from a single seed. """ -default_solver() = DPWSolver(; n_iterations = 100_000, tree_in_info = true) +default_solver(rng::AbstractRNG = default_rng()) = + DPWSolver(; n_iterations = 100_000, tree_in_info = true, rng) # minimize the expected experimental cost while ensuring the uncertainty remains below a specified threshold. include("UncertaintyReductionMDP.jl") diff --git a/src/GenerativeDesigns/UncertaintyReductionMDP.jl b/src/GenerativeDesigns/UncertaintyReductionMDP.jl index 077213f..5adf60d 100644 --- a/src/GenerativeDesigns/UncertaintyReductionMDP.jl +++ b/src/GenerativeDesigns/UncertaintyReductionMDP.jl @@ -5,6 +5,14 @@ Structure that parametrizes the experimental decision-making process. It is used In this experimental setup, our objective is to minimize the expected experimental cost while ensuring the uncertainty remains below a specified threshold. +!!! note "Relationship to the paper" + The uncertainty measure returned by [`Variance`](@ref) / [`Entropy`](@ref) is + *normalized* — it is the fraction of the prior uncertainty — so `threshold` is + interpreted on `[0, 1]` (this is why `efficient_designs` sweeps thresholds over + `range(0, 1, …)`). The paper's terminal information term is realized here as a + **hard terminal condition** (`uncertainty ≤ threshold`) rather than as an additive + reward penalty; the objective is otherwise equivalent. + Internally, a state of the decision process is modeled as a tuple `(evidence::Evidence, [total accumulated monetary cost, total accumulated execution time])`. The per-step reward is the negative marginal cost incurred at that step, i.e. @@ -79,8 +87,16 @@ struct UncertaintyReductionMDP{S, U} <: POMDPs.MDP{State, Vector{String}} state = State((evidence, Tuple(zeros(2)))) # check if `sampler`, `uncertainty` are compatible - @assert hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) """`sampler` must implement a method accepting `(evidence, readout features, rng)` as its arguments.""" - @assert hasmethod(uncertainty, Tuple{Evidence}) """`uncertainty` must implement a method accepting `evidence` as its argument.""" + hasmethod(sampler, Tuple{Evidence, Vector{String}, AbstractRNG}) || throw( + ArgumentError( + "`sampler` must implement a method accepting `(evidence, readout features, rng)` as its arguments.", + ), + ) + hasmethod(uncertainty, Tuple{Evidence}) || throw( + ArgumentError( + "`uncertainty` must implement a method accepting `evidence` as its argument.", + ), + ) # actions and their costs costs = Dict{String, ActionCost}( @@ -129,7 +145,15 @@ function POMDPs.actions(m::UncertaintyReductionMDP, state) return !isempty(feats) && !any(f -> haskey(state.evidence, f), feats) end - return if !isempty(all_actions) && (length(state.evidence) < m.max_experiments) + # Count *completed experiments* (all of an experiment's features present), not + # raw evidence entries, so prior evidence / multi-feature experiments don't + # mis-count against `max_experiments`. + n_completed = count(keys(m.costs)) do a + feats = m.costs[a].features + !isempty(feats) && all(f -> haskey(state.evidence, f), feats) + end + + return if !isempty(all_actions) && (n_completed < m.max_experiments) collect(powerset(all_actions, 1, m.max_parallel)) else [[eox]] @@ -228,7 +252,7 @@ function efficient_design( uncertainty, threshold, evidence = Evidence(), - solver = default_solver(), + solver = nothing, repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -254,12 +278,12 @@ function efficient_design( (; monetary_cost = 0.0, time = 0.0), ) else - # planner - planner = solve(solver, mdp) + # planner (fresh solver seeded by `rng` unless the caller supplied one) + planner = solve(isnothing(solver) ? default_solver(rng) : solver, mdp) action, info = action_info(planner, mdp.initial_state) if repetitions > 0 - queue = [Sim(mdp, planner; rng) for _ in 1:repetitions] + queue = [Sim(mdp, planner; rng = Xoshiro(rand(rng, UInt64))) for _ in 1:repetitions] stats = run_parallel(queue) do _, hist monetary_cost, time = hist[end][:s].costs @@ -362,7 +386,7 @@ function efficient_designs( uncertainty, thresholds, evidence = Evidence(), - solver = default_solver(), + solver = nothing, repetitions = 0, realized_uncertainty = false, mdp_options = (;), @@ -372,6 +396,10 @@ function efficient_designs( designs = [] for threshold in range(0.0, 1.0, thresholds) @info "Current threshold level : $threshold" + # Each threshold gets an independent rng (derived from the master `rng`) + # and, unless the caller supplied a solver, a fresh solver seeded by it. + inner_rng = Xoshiro(rand(rng, UInt64)) + inner_solver = isnothing(solver) ? default_solver(inner_rng) : solver push!( designs, efficient_design( @@ -380,11 +408,11 @@ function efficient_designs( uncertainty, threshold, evidence, - solver, + solver = inner_solver, repetitions, realized_uncertainty, mdp_options, - rng, + rng = inner_rng, ), ) end diff --git a/src/GenerativeDesigns/distancebased.jl b/src/GenerativeDesigns/distancebased.jl index 5f4077e..f491508 100644 --- a/src/GenerativeDesigns/distancebased.jl +++ b/src/GenerativeDesigns/distancebased.jl @@ -6,11 +6,18 @@ This returns an anonymous function `(x, col; prior) -> λ * (x .- col).^2 / σ`. If `standardize` is set to `true`, `σ` represents `col`'s variance calculated in relation to `prior`, otherwise `σ` equals one. """ function QuadraticDistance(; λ = 1, standardize = true) - σ = nothing - + # NOTE: `σ` is recomputed on every call rather than memoized. A memoized `σ` + # would be reused across columns if the same closure instance is shared + # (e.g. `Dict(c => qd for c in cols)` with a single `qd`), applying the + # first column's variance to the others. Recomputing keeps each call correct. return function (x, col; prior = ones(length(col))) - if isnothing(σ) - σ = standardize ? var(col, Weights(prior); corrected = false) : 1 + σ = standardize ? var(col, Weights(prior); corrected = false) : 1 + if standardize && σ <= 0 + throw( + ArgumentError( + "QuadraticDistance: column has zero variance under the prior (constant column); cannot standardize. Pass `standardize = false` or remove the constant column.", + ), + ) end return λ * (x .- col) .^ 2 / σ @@ -20,17 +27,23 @@ end """ DiscreteDistance(; λ=1) -Return an anonymous function `(x, col) -> λ * (x .== col)`. +Return an anonymous function `(x, col) -> λ * (x .!= col)`. + +This is a proper distance: matching entries contribute `0` (nearest) and +mismatching entries contribute `λ`, consistent with [`QuadraticDistance`](@ref). +The collated distances are passed through the `similarity` functional +(e.g. [`Exponential`](@ref)), so rows matching the evidence receive the +highest similarity weight. """ DiscreteDistance(; λ = 1) = function (x, col; _...) - return map(y -> y == x ? λ : 0.0, col) + return map(y -> y == x ? 0.0 : λ, col) end # Default similarity functional """ - Exponential(; λ=1) + Exponential(; λ=1/2) -Return an anonymous function `x -> exp(-λ * sum(x; init=0))`. +Return an anonymous function `x -> exp(-λ * x)`. """ Exponential(; λ = 1 / 2) = x -> exp(-λ * x) @@ -72,7 +85,11 @@ it returns an internal function of `weights` that computes the fraction of infor """ function Entropy() return function (labels; prior) - @assert elscitype(labels) <: Multiclass "labels must be of `Multiclass` scitype, but `elscitype(labels)=$(elscitype(labels))`" + elscitype(labels) <: Multiclass || throw( + ArgumentError( + "labels must be of `Multiclass` scitype, but `elscitype(labels)=$(elscitype(labels))`", + ), + ) initial = compute_entropy(labels; weights = prior) initial > 0 || throw( ArgumentError( @@ -232,14 +249,14 @@ A named tuple with the following fields: - `sampler`: a function of `(evidence, features, rng)`, in which `evidence` denotes the current experimental evidence, `features` represent the set of features we want to sample from, and `rng` is a random number generator; it returns a dictionary mapping the features to outcomes. - `uncertainty`: a function of `evidence`; it returns the measure of variance or uncertainty about the target variable, conditioned on the experimental evidence acquired so far. - - `weights`: a function of `evidence`; it returns probabilities (posterior) acrss the rows in `data`. + - `weights`: a function of `evidence`; it returns probabilities (posterior) across the rows in `data`. # Arguments - `data`: a dataframe with historical data. - `target`: target column name or a vector of target columns names. -# Keyword Argumets +# Keyword Arguments - `uncertainty`: a function that takes the subdataframe containing columns in targets along with prior, and returns an anonymous function taking a single argument (a probability vector over observations) and returns an uncertainty measure over targets. - `similarity`: a function that, for each row, takes distances between `row[col]` and `readout[col]`, and returns a non-negative probability mass for the row. @@ -312,11 +329,13 @@ function DistanceBased( # If "importance weight" is a function, apply it to the column to get a numeric vector. for (colname, val) in importance_weights - push!( - weights, - string(colname) => - (val isa Function ? map(x -> val(x), data[!, colname]) : val), + wv = val isa Function ? map(x -> val(x), data[!, colname]) : val + length(wv) == nrow(data) || throw( + ArgumentError( + "`importance_weights` for column `$(colname)` has length $(length(wv)), but `data` has $(nrow(data)) rows.", + ), ) + push!(weights, string(colname) => wv) end # Convert "desirable ranges" into importance weights. diff --git a/src/StaticDesigns/StaticDesigns.jl b/src/StaticDesigns/StaticDesigns.jl index 8e8f9d7..bd3ee1a 100644 --- a/src/StaticDesigns/StaticDesigns.jl +++ b/src/StaticDesigns/StaticDesigns.jl @@ -6,6 +6,7 @@ using Combinatorics: powerset using POMDPs using POMDPTools: Deterministic using MCTS +using Random using ..CEEDesigns: front @@ -86,7 +87,7 @@ function evaluate_experiments( foreach( x -> append!( features, - experiments[x] isa Pair ? experiments[x][2] : experiments[x], + string.(experiments[x] isa Pair ? experiments[x][2] : experiments[x]), ), exp_set, ) @@ -159,7 +160,7 @@ function evaluate_experiments( foreach( x -> append!( features, - experiments[x] isa Pair ? experiments[x][2] : experiments[x], + string.(experiments[x] isa Pair ? experiments[x][2] : experiments[x]), ), exp_set, ) @@ -219,6 +220,7 @@ function efficient_designs( max_parallel::Int = 1, tradeoff = (1, 0), mdp_kwargs = default_mdp_kwargs, + rng::AbstractRNG = Random.default_rng(), ) where {T, S} experimental_costs = Dict(e => v isa Pair ? v[1] : v for (e, v) in experiments) @@ -238,7 +240,13 @@ function efficient_designs( # lock to prevent race condition lk = ReentrantLock() - Threads.@threads for design in collect(evals) + tasks = collect(evals) + # Pre-seed one independent rng per task so the threaded arrangement search is + # reproducible and free of GLOBAL_RNG contention across worker threads. + seeds = rand(rng, UInt64, length(tasks)) + + Threads.@threads for i in eachindex(tasks) + design = tasks[i] arrangement = optimal_arrangement( experimental_costs, evals, @@ -246,6 +254,7 @@ function efficient_designs( max_parallel, tradeoff, mdp_kwargs, + rng = Xoshiro(seeds[i]), ) lock(lk) do diff --git a/src/StaticDesigns/arrangements.jl b/src/StaticDesigns/arrangements.jl index 9d94e15..7f74ec3 100644 --- a/src/StaticDesigns/arrangements.jl +++ b/src/StaticDesigns/arrangements.jl @@ -69,6 +69,7 @@ function optimal_arrangement( max_parallel = 1, tradeoff = (1, 0), mdp_kwargs = default_mdp_kwargs, + rng::AbstractRNG = Random.default_rng(), ) where {T} # validate that `evals` covers every non-terminal state the MDP rollout may visit. # the rollout visits all subsets of `experiments` of size 0..|experiments|-1, so each @@ -107,7 +108,7 @@ function optimal_arrangement( mdp = ArrangementMDP(; experiments, experimental_costs, evals, max_parallel, tradeoff) - solver = DPWSolver(; mdp_kwargs...) + solver = DPWSolver(; mdp_kwargs..., rng) planner = solve(solver, mdp) monetary_cost = time = 0.0 diff --git a/src/ensemble_fronts.jl b/src/ensemble_fronts.jl index b10cb89..f66c2f5 100644 --- a/src/ensemble_fronts.jl +++ b/src/ensemble_fronts.jl @@ -5,11 +5,17 @@ using DataFrames Flatten a vector of Pareto fronts (one per ensemble run) into `df::DataFrame`: DataFrame with columns Threshold, Action_Set, -Frequency, Average_Utility +Frequency, Average_Cost """ function ensemble_to_dataframe(runs) counts = Dict{Tuple{Float64, String}, Tuple{Int, Float64}}() for front in runs + # De-duplicate identical (threshold, actions) points *within a single run* + # before counting. With `realized_uncertainty = true`, several + # already-terminal thresholds can collapse onto identical points that + # `front()` (correctly) keeps; counting each would inflate this run's + # frequency (and the MLASP vote). Each run contributes at most +1 per key. + seen = Dict{Tuple{Float64, String}, Float64}() for design in front threshold = design[1][2] actions = if haskey(design[2], :arrangement) @@ -22,17 +28,20 @@ function ensemble_to_dataframe(runs) else "" end - utility = design[1][1] + cost = design[1][1] key = (threshold, actions) + haskey(seen, key) || (seen[key] = cost) + end + for (key, cost) in seen prev = get(counts, key, (0, 0.0)) - counts[key] = (prev[1] + 1, prev[2] + utility) + counts[key] = (prev[1] + 1, prev[2] + cost) end end rows = [ ( Threshold = k[1], Action_Set = k[2], - Average_Utility = v[2] / v[1], + Average_Cost = v[2] / v[1], Frequency = v[1], ) for (k, v) in counts ] @@ -58,7 +67,7 @@ function _compute_dispersion(df_copy::DataFrame) dispersion_offsets = Dict{Float64, Vector{Float64}}() action_counts = Dict{Float64, Int}() - x_range = maximum(df_copy.Average_Utility) - minimum(df_copy.Average_Utility) + x_range = maximum(df_copy.Average_Cost) - minimum(df_copy.Average_Cost) for group in threshold_groups threshold = group.Threshold[1] @@ -95,7 +104,7 @@ function _compute_dispersion(df_copy::DataFrame) # threshold while iterating `eachrow(df_copy)` corresponds to position n in # `dispersion_offsets[threshold]`, because `groupby` preserves within-group # row order. An explicit per-threshold counter avoids float `==` lookups - # that would collapse duplicate `Average_Utility` rows onto the same offset + # that would collapse duplicate `Average_Cost` rows onto the same offset # slot. x_coords = Float64[] threshold_indices = Dict{Float64, Int}() @@ -113,9 +122,9 @@ function _compute_dispersion(df_copy::DataFrame) if haskey(dispersion_offsets, threshold) && point_index <= length(dispersion_offsets[threshold]) offset = dispersion_offsets[threshold][point_index] - push!(x_coords, row.Average_Utility + offset) + push!(x_coords, row.Average_Cost + offset) else - push!(x_coords, row.Average_Utility) + push!(x_coords, row.Average_Cost) end end @@ -129,7 +138,7 @@ Create Pareto front visualization for ensemble results with viridis color scheme # Arguments - - `df::DataFrame`: DataFrame with columns Threshold, Action_Set, Frequency, Average_Utility + - `df::DataFrame`: DataFrame with columns Threshold, Action_Set, Frequency, Average_Cost - `tau::Float64`: Belief threshold value for labeling - `show_annotations::Bool`: Whether to show action set annotations (default: true) @@ -162,7 +171,7 @@ function plot_ensemble_pareto( normalization::Symbol = :global, ) # Guard against empty input. Several downstream calls (e.g. - # `maximum(values(threshold_counts))`, `maximum(df_copy.Average_Utility)`) + # `maximum(values(threshold_counts))`, `maximum(df_copy.Average_Cost)`) # would otherwise raise an opaque `ArgumentError` from `maximum`/`argmax`. isempty(df) && throw( ArgumentError( @@ -494,7 +503,7 @@ function plot_ensemble_pareto( push!( annotations_to_place, ( - x = row.Average_Utility / scale_factor, # Scale x-coordinate + x = row.Average_Cost / scale_factor, # Scale x-coordinate y = row.Threshold, label = label_text, ), diff --git a/src/fronts.jl b/src/fronts.jl index 7638a17..5d18486 100644 --- a/src/fronts.jl +++ b/src/fronts.jl @@ -118,7 +118,7 @@ function front( end """ - plot_front(designs; grad=cgrad(:Paired_12), xlabel, ylabel, labels=get_labels(designs)) + plot_front(designs; grad=cgrad(:Paired_12), xlabel, ylabel, labels=make_labels(designs)) Render scatter plot of efficient designs, as returned from `efficient_designs`. @@ -139,6 +139,9 @@ function plot_front( ylabel = "information measure", labels = make_labels(designs), ) + isempty(designs) && throw( + ArgumentError("plot_front requires at least one design (got an empty collection)."), + ) xs = map(x -> x[1][1], designs) ys = map(x -> x[1][2], designs) diff --git a/test/GenerativeDesigns/test.jl b/test/GenerativeDesigns/test.jl index 1216537..0a78f73 100644 --- a/test/GenerativeDesigns/test.jl +++ b/test/GenerativeDesigns/test.jl @@ -1,6 +1,9 @@ # test with the sum distances include("test_distances_sum.jl") +# regression tests for pre-JOSS review fixes (S1/C1/C2/M1/M4/N3) +include("test_review_fixes.jl") + # test with the squared Mahalanobis distance include("test_mahalanobis.jl") diff --git a/test/GenerativeDesigns/test_efficient_value.jl b/test/GenerativeDesigns/test_efficient_value.jl index 43ed3e7..fe8a41f 100644 --- a/test/GenerativeDesigns/test_efficient_value.jl +++ b/test/GenerativeDesigns/test_efficient_value.jl @@ -80,7 +80,7 @@ end # (2) A `value` with the wrong (Evidence, Vector{Float64}) signature must be # rejected. The buggy assertion would have accepted this even though the # MDP itself would later pass `state.costs::NTuple{2, Float64}` and crash. -@test_throws AssertionError EfficientValueMDP( +@test_throws ArgumentError EfficientValueMDP( experiments; sampler, value = value_vector, diff --git a/test/GenerativeDesigns/test_review_fixes.jl b/test/GenerativeDesigns/test_review_fixes.jl new file mode 100644 index 0000000..2c99175 --- /dev/null +++ b/test/GenerativeDesigns/test_review_fixes.jl @@ -0,0 +1,116 @@ +using Test +using DataFrames +using ScientificTypes +import POMDPs +using CEEDesigns, CEEDesigns.GenerativeDesigns + +@testset "Pre-JOSS review fixes" begin + + # ---- S1: DiscreteDistance is a proper distance (match -> 0, mismatch -> λ) ---- + @testset "S1 DiscreteDistance direction" begin + d = DiscreteDistance(; λ = 1) + @test d("M", ["M", "F", "M"]) == [0.0, 1.0, 0.0] + + data = DataFrame(; target = [1.0, 2.0, 3.0, 4.0], Sex = ["M", "M", "F", "F"]) + coerce!(data, :Sex => Multiclass) + (; weights) = DistanceBased( + data; + target = "target", + uncertainty = Variance(), + similarity = Exponential(; λ = 1.0), + distance = Dict("Sex" => DiscreteDistance()), + ) + w = weights(Evidence("Sex" => "M")) + # rows 1,2 match the query (Sex=M) -> must be weighted strictly higher + @test sum(w[1:2]) > sum(w[3:4]) + end + + # ---- C1: terminal_condition accepts Symbol keys (no runtime crash) ---- + @testset "C1 terminal_condition Symbol keys" begin + data = DataFrame(; target = [0.2, 0.7, 0.9], x = [1, 2, 3]) + weights_fn = ev -> [0.2, 0.3, 0.5] + unc = ev -> 0.0 + sampler = (ev, f, rng) -> Dict(f[1] => 1.0) + costs = Dict("A" => 1.0) + mdp = ConditionalUncertaintyReductionMDP( + costs; + sampler, + uncertainty = unc, + threshold = 0.1, + weights = weights_fn, + data = data, + terminal_condition = (Dict(:target => [0.5, 1.0]), 0.75), # SYMBOL key + ) + # admissible mass = 0.3 + 0.5 = 0.8 >= 0.75 -> terminal, and must not throw + @test POMDPs.isterminal(mdp, mdp.initial_state) == true + end + + # ---- C2: Union{Missing,Real} constraint column accepted; missing excluded ---- + @testset "C2 Union{Missing} constraint column" begin + data = DataFrame(; target = Union{Missing, Float64}[0.2, 0.7, missing], x = [1, 2, 3]) + weights_fn = ev -> [0.2, 0.3, 0.5] + unc = ev -> 0.0 + sampler = (ev, f, rng) -> Dict(f[1] => 1.0) + costs = Dict("A" => 1.0) + mdp = ConditionalUncertaintyReductionMDP( + costs; + sampler, + uncertainty = unc, + threshold = 0.1, + weights = weights_fn, + data = data, + terminal_condition = (Dict("target" => [0.5, 1.0]), 0.2), + ) + @test mdp isa ConditionalUncertaintyReductionMDP + # row1=0.2 (out), row2=0.7 (in -> 0.3), row3=missing (excluded) + p = conditional_likelihood( + Evidence(); + compute_weights = weights_fn, + hist_data = data, + target_condition = Dict("target" => (0.5, 1.0)), + ) + @test p ≈ 0.3 + end + + # ---- M1: max_experiments counts experiments, not evidence entries ---- + @testset "M1 max_experiments counts experiments" begin + sampler = (ev, f, rng) -> Dict(f[1] => 1.0) + unc = ev -> 1.0 # never terminal via uncertainty + costs = Dict("ECG" => ((1.0, 0.5) => ["MaxHR", "RestingECG"])) + mdp = UncertaintyReductionMDP( + costs; + sampler, + uncertainty = unc, + threshold = 0.5, + evidence = Evidence("Age" => 35, "Sex" => "M"), + max_experiments = 2, + ) + acts = POMDPs.actions(mdp, mdp.initial_state) + @test acts != [["EOX"]] # 2 prior-evidence entries must NOT exhaust budget + @test any(a -> "ECG" in a, acts) + end + + # ---- M4: constant column raises a clear ArgumentError ---- + @testset "M4 constant column ArgumentError" begin + data = DataFrame(; target = [1.0, 2.0, 3.0], c = [5.0, 5.0, 5.0]) + (; weights) = DistanceBased( + data; + target = "target", + uncertainty = Variance(), + distance = Dict("c" => QuadraticDistance()), + ) + @test_throws ArgumentError weights(Evidence("c" => 5.0)) + end + + # ---- N3: importance_weights length is validated ---- + @testset "N3 importance_weights length" begin + data = DataFrame(; target = [1.0, 2.0, 3.0], c = [5.0, 6.0, 7.0]) + @test_throws ArgumentError DistanceBased( + data; + target = "target", + uncertainty = Variance(), + importance_weights = Dict("c" => [1.0, 2.0]), # length 2 != 3 rows + ) + end + +end diff --git a/test/StaticDesigns/test.jl b/test/StaticDesigns/test.jl index ce9f3cd..7a4c70b 100644 --- a/test/StaticDesigns/test.jl +++ b/test/StaticDesigns/test.jl @@ -1,4 +1,5 @@ using Random: seed! +using Random: Xoshiro using CEEDesigns, CEEDesigns.StaticDesigns using CSV, DataFrames @@ -144,3 +145,28 @@ truncated_evals = Dict( Set(keys(experiments_simple)); mdp_kwargs = (; n_iterations = 50, depth = 5, exploration_constant = 3.0), ) + +## S3 regression: seeded rng => reproducible arrangement (no GLOBAL_RNG contention). +let + repro_costs = Dict("A" => 1.0, "B" => 2.0, "C" => 3.0) + repro_evals = + Dict{Set{String}, NamedTuple{(:loss, :filtration), Tuple{Float64, Float64}}}( + Set{String}() => (; loss = 1.0, filtration = 1.0), + Set(["A"]) => (; loss = 0.7, filtration = 0.9), + Set(["B"]) => (; loss = 0.6, filtration = 0.8), + Set(["C"]) => (; loss = 0.5, filtration = 0.7), + Set(["A", "B"]) => (; loss = 0.4, filtration = 0.6), + Set(["A", "C"]) => (; loss = 0.4, filtration = 0.6), + Set(["B", "C"]) => (; loss = 0.4, filtration = 0.6), + ) + full = Set(["A", "B", "C"]) + mdp_kwargs = (; n_iterations = 500, depth = 5, exploration_constant = 3.0) + a1 = CEEDesigns.StaticDesigns.optimal_arrangement( + repro_costs, repro_evals, full; mdp_kwargs, rng = Xoshiro(123), + ) + a2 = CEEDesigns.StaticDesigns.optimal_arrangement( + repro_costs, repro_evals, full; mdp_kwargs, rng = Xoshiro(123), + ) + @test a1.arrangement == a2.arrangement + @test a1.combined_cost == a2.combined_cost +end diff --git a/test/test_ensemble_fronts.jl b/test/test_ensemble_fronts.jl index 799758e..202e0d6 100644 --- a/test/test_ensemble_fronts.jl +++ b/test/test_ensemble_fronts.jl @@ -38,16 +38,32 @@ end @test nrow(df2) == 1 @test df2.Action_Set[1] == "A,B" @test df2.Frequency[1] == 2 - @test df2.Average_Utility[1] == 1.5 + @test df2.Average_Cost[1] == 1.5 end -@testset "_compute_dispersion distributes duplicate Average_Utility (M11)" begin - # Two rows with identical Average_Utility under the same threshold must NOT +@testset "ensemble_to_dataframe de-duplicates within a run (M2)" begin + # Two identical collapsed points in ONE run (same threshold, same empty + # action set) must contribute frequency 1, not 2. + dup_design = ((0.0, 0.25), (; monetary_cost = 0.0, time = 0.0)) # no :arrangement => "" + df = ensemble_to_dataframe([[dup_design, dup_design]]) + @test nrow(df) == 1 + @test df.Frequency[1] == 1 + @test df.Threshold[1] == 0.25 + @test df.Action_Set[1] == "" + + # The same key across TWO separate runs still counts as frequency 2. + df2 = ensemble_to_dataframe([[dup_design], [dup_design]]) + @test nrow(df2) == 1 + @test df2.Frequency[1] == 2 +end + +@testset "_compute_dispersion distributes duplicate Average_Cost (M11)" begin + # Two rows with identical Average_Cost under the same threshold must NOT # both be assigned the slot-1 offset. df = DataFrame(; Threshold = [0.5, 0.5, 0.5], Action_Set = ["A", "B", "C"], - Average_Utility = [1.0, 1.0, 2.0], + Average_Cost = [1.0, 1.0, 2.0], Frequency = [3, 1, 2], ) @@ -59,7 +75,7 @@ end # `range(-max_offset, max_offset; length = n_actions)` of length 3). @test length(unique(dispersion_offsets[0.5])) == 3 # The x_coords vector must produce three distinct values, even though two - # rows share Average_Utility. This is the regression: the old findfirst + + # rows share Average_Cost. This is the regression: the old findfirst + # float-equality logic collapsed both duplicates onto offset slot 1, so # x_coords[1] == x_coords[2]. With the per-threshold counter fix they # consume slots 1 and 2 respectively. @@ -71,7 +87,7 @@ end empty_df = DataFrame(; Threshold = Float64[], Action_Set = String[], - Average_Utility = Float64[], + Average_Cost = Float64[], Frequency = Int[], ) @test_throws ArgumentError plot_ensemble_pareto(empty_df, 0.5)