From 6b6b8e854f4275079881feffd921293f353f8dbd Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 20 Jun 2026 06:27:41 -0400 Subject: [PATCH 1/5] Fix master test reds: BBO loss return, ensemble solve API, FlexiChains chain access Three independent breakages from upstream SciML updates were turning the Core and Datafit CI groups red on master: 1. Optimization loss return type (Datafit group). `l2loss`/`relative_l2loss` returned `(tot_loss, sol)`. OptimizationBBO's BBO wrapper now strictly requires the objective to return a `Float64`, so `global_datafit` errored with "fitness function does NOT return the expected fitness type Float64". The `sol` element of the tuple was never consumed by any caller, so the loss functions now return only the scalar `tot_loss`. NLopt's `datafit` path (which tolerated the tuple by taking `first`) is unaffected. 2. EnsembleProblem / ensemble solve (Core group). The vector-of-problems form `EnsembleProblem([prob1, prob2, prob3])` is deprecated and broken in current SciMLBase (the default prob_func passes the whole vector to the per-trajectory solve), so `solve(enprob; saveat=1)` failed with a MethodError on `init(::Vector{ODEProblem}, ::Nothing)`. Migrated `bayesian_ensemble` and the ensemble test to the modern `prob_func = (prob, ctx) -> probs[ctx.sim_id]` form, pass an explicit algorithm (`Tsit5()`) and `trajectories`, and access per-trajectory solutions via `sol.u[i]` (EnsembleSolution's `length`/symbolic `getindex` now flatten across all timepoints). `ensemble_weights` updated to use `sol.u` accordingly. 3. Turing chain backend (Core + Datafit groups). Turing now returns a `FlexiChains.VNChain`, which no longer supports `chain["pprior[i]"]` string indexing. `bayesian_datafit` now extracts posterior samples with `chain[@varname(pprior[i])]`. Verified locally on Julia 1.12 against the master dependency set: - global_datafit (BBO) recovers [2/3, 4/3, 1, 1]; datafit (NLopt) still passes. - prob_func ensemble solve produces a Vector{ODESolution}; ensemble_weights runs. - bayesian_datafit returns per-parameter posterior sample vectors with reduced variance vs the prior (the Datafit test's assertion). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- src/datafit.jl | 12 ++++++------ src/ensemble.jl | 6 ++++-- test/ensemble.jl | 19 +++++++++++-------- 3 files changed, 21 insertions(+), 16 deletions(-) diff --git a/src/datafit.jl b/src/datafit.jl index b40fb0f..ab6a3c5 100644 --- a/src/datafit.jl +++ b/src/datafit.jl @@ -6,7 +6,7 @@ function l2loss(pvals, (prob, pkeys, t, data)::Tuple{Vararg{Any, 4}}) for pairs in data tot_loss += sum((sol[pairs.first] .- pairs.second) .^ 2) end - return tot_loss, sol + return tot_loss end function l2loss(pvals, (prob, pkeys, data)::Tuple{Vararg{Any, 3}}) @@ -22,7 +22,7 @@ function l2loss(pvals, (prob, pkeys, data)::Tuple{Vararg{Any, 3}}) for i in 1:length(ts) tot_loss += sum((sol(ts[i]; idxs = datakeys[i]) .- timeseries[i]) .^ 2) end - return tot_loss, sol + return tot_loss end function relative_l2loss(pvals, (prob, pkeys, t, data)::Tuple{Vararg{Any, 4}}) @@ -33,7 +33,7 @@ function relative_l2loss(pvals, (prob, pkeys, t, data)::Tuple{Vararg{Any, 4}}) for pairs in data tot_loss += sum(((sol[pairs.first] .- pairs.second) ./ sol[pairs.first]) .^ 2) end - return tot_loss, sol + return tot_loss end function relative_l2loss(pvals, (prob, pkeys, data)::Tuple{Vararg{Any, 3}}) @@ -50,7 +50,7 @@ function relative_l2loss(pvals, (prob, pkeys, data)::Tuple{Vararg{Any, 3}}) vals = sol(ts[i]; idxs = datakeys[i]) tot_loss += sum(((vals .- timeseries[i]) ./ vals) .^ 2) end - return tot_loss, sol + return tot_loss end """ @@ -308,7 +308,7 @@ function bayesian_datafit( progress = true ) return [ - Pair(p[i].first, collect(chain["pprior[" * string(i) * "]"])[:]) + Pair(p[i].first, collect(chain[@varname(pprior[i])])[:]) for i in eachindex(p) ] end @@ -333,7 +333,7 @@ function bayesian_datafit( progress = true ) return [ - Pair(p[i].first, collect(chain["pprior[" * string(i) * "]"])[:]) + Pair(p[i].first, collect(chain[@varname(pprior[i])])[:]) for i in eachindex(p) ] end diff --git a/src/ensemble.jl b/src/ensemble.jl index ded0389..a4480d0 100644 --- a/src/ensemble.jl +++ b/src/ensemble.jl @@ -19,7 +19,7 @@ dataset on which the ensembler should be trained on. function ensemble_weights(sol::EnsembleSolution, data_ensem) obs = first.(data_ensem) predictions = reduce( - vcat, reduce(hcat, [sol[i][s] for i in 1:length(sol)]) for s in obs + vcat, reduce(hcat, [sol.u[i][s] for i in 1:length(sol.u)]) for s in obs ) data = reduce( vcat, @@ -56,5 +56,7 @@ function bayesian_ensemble( @info "$(length(all_probs)) total models" - return enprob = EnsembleProblem(all_probs) + return enprob = EnsembleProblem( + all_probs[1]; prob_func = (prob, ctx) -> all_probs[ctx.sim_id] + ) end diff --git a/test/ensemble.jl b/test/ensemble.jl index 08bfa7e..39f9da3 100644 --- a/test/ensemble.jl +++ b/test/ensemble.jl @@ -52,15 +52,16 @@ eqs = [ @mtkbuild sys3 = ODESystem(eqs, t) prob3 = ODEProblem(sys3, [], tspan); -enprob = EnsembleProblem([prob, prob2, prob3]) +probs = [prob, prob2, prob3] +enprob = EnsembleProblem(probs[1]; prob_func = (prob, ctx) -> probs[ctx.sim_id]) -sol = solve(enprob; saveat = 1); +sol = solve(enprob, Tsit5(); saveat = 1, trajectories = length(probs)); weights = [0.2, 0.5, 0.3] -fullS = vec(sum(stack(weights .* sol[S, :]), dims = 2)) -fullI = vec(sum(stack(weights .* sol[I, :]), dims = 2)) -fullR = vec(sum(stack(weights .* sol[R, :]), dims = 2)) +fullS = vec(sum(stack(weights .* [sol.u[i][S] for i in 1:length(sol.u)]), dims = 2)) +fullI = vec(sum(stack(weights .* [sol.u[i][I] for i in 1:length(sol.u)]), dims = 2)) +fullR = vec(sum(stack(weights .* [sol.u[i][R] for i in 1:length(sol.u)]), dims = 2)) t_train = 0:14 data_train = [ @@ -81,14 +82,16 @@ data_forecast = [ R => (t_forecast, fullR), ] -sol = solve(enprob; saveat = t_ensem); +sol = solve(enprob, Tsit5(); saveat = t_ensem, trajectories = length(probs)); @test ensemble_weights(sol, data_ensem) ≈ [0.2, 0.5, 0.3] -probs = [prob, prob2, prob3] ps = [[β => Uniform(0.01, 10.0), γ => Uniform(0.01, 10.0)] for i in 1:3] datas = [data_train, data_train, data_train] enprobs = bayesian_ensemble(probs, ps, datas) -sol = solve(enprobs; saveat = t_ensem); +sol = solve( + enprobs, Tsit5(); saveat = t_ensem, + trajectories = length(enprobs.prob_func.all_probs) +); ensemble_weights(sol, data_ensem) From 1b8a9ffe62ca9ec97c30851afb06dcca4c24e301 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 20 Jun 2026 08:42:19 -0400 Subject: [PATCH 2/5] Fix Core/Datafit/Downgrade reds: DynamicPPL @addlogprob!, cross-version ensemble prob_func Builds on the prior fixes in this branch (BBO scalar loss return, FlexiChains @varname extraction). Three remaining breakages from upstream SciML/Turing churn: 1. Datafit group (DynamicPPL 0.41). bayesianODE rejected non-successful solves with `Turing.DynamicPPL.acclogp!!(__varinfo__, -Inf)`, which no longer has a method when __varinfo__ is the AD-time `OnlyAccsVarInfo{...Dual...}`: MethodError: no method matching acclogp!!(::OnlyAccsVarInfo{...}, ::Float64) Replaced with the canonical sample-rejection idiom `@addlogprob! (; loglikelihood = -Inf)`, which dispatches correctly on the LogLikelihood accumulator under ForwardDiff. 2/3. Core + Downgrade groups (EnsembleProblem prob_func arity). The prior `(prob, ctx) -> probs[ctx.sim_id]` form only works on SciMLBase 3.x (2-arg `prob_func(prob, ctx)`). On the Downgrade floor (SciMLBase 2.55) EnsembleProblem still calls the legacy 3-arg `prob_func(prob, i, repeat)`, so it errored with MethodError: no method matching (::#1#2)(::ODEProblem, ::Int64, ::Int64) Introduced `EnsembleProbForwarder(all_probs)`, a callable supporting both interfaces (integer index and `ctx.sim_id`), used as `bayesian_ensemble`'s prob_func. Storing `all_probs` lets callers read the trajectory count via `enprob.prob_func.all_probs` (the access the ensemble test already uses). The ensemble test's inline prob_func is likewise made arity-robust. Verified locally on Julia 1.12 against the master dependency set (SciMLBase 3.21, DynamicPPL 0.41.8, Turing 0.45, OptimizationBBO 0.4.7): - simple EnsembleProblem solve(enprob, Tsit5(); trajectories) recovers weights [0.2, 0.5, 0.3]; EnsembleProbForwarder dispatches on both (prob, i, repeat) and (prob, ctx); ensemble_weights runs. - bayesian_datafit (both t and (t, timeseries) forms) returns per-parameter posterior sample vectors with variance reduced vs the prior (the Datafit assertion). Not fixed here (upstream, reported separately): QA group's `Aqua.find_persistent_tasks_deps` fails ONLY on Julia 1.12 (passes on lts) because Pkg now honors OptimizationBBO 0.4.4-0.4.7's released relative-path `[sources] OptimizationBase = {path = "../OptimizationBase"}` and errors with "expected package OptimizationBase to exist at path .../OptimizationBBO/OptimizationBase". Same OptimizationBBO version passes on Julia 1.10. Capping OptimizationBBO < 0.4.4 would drag SciMLBase 3->2, ModelingToolkit 11->9, OptimizationBase 5->2 (major regression), so the correct fix is upstream (stop shipping [sources] in releases). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- src/datafit.jl | 4 ++-- src/ensemble.jl | 19 ++++++++++++++++++- test/ensemble.jl | 4 +++- 3 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/datafit.jl b/src/datafit.jl index ab6a3c5..8b7208c 100644 --- a/src/datafit.jl +++ b/src/datafit.jl @@ -223,7 +223,7 @@ Turing.@model function bayesianODE(prob, t, pdist, pkeys, data, noise_prior) prob = remake(prob, tspan = (prob.tspan[1], t[end]), p = Pair.(pkeys, pprior)) sol = solve(prob, saveat = t) if !SciMLBase.successful_retcode(sol) - Turing.DynamicPPL.acclogp!!(__varinfo__, -Inf) + Turing.@addlogprob! (; loglikelihood = -Inf) return nothing end for i in eachindex(data) @@ -249,7 +249,7 @@ Turing.@model function bayesianODE( prob = remake(prob, tspan = (prob.tspan[1], lastt), p = Pair.(pkeys, pprior)) sol = solve(prob) if !SciMLBase.successful_retcode(sol) - Turing.DynamicPPL.acclogp!!(__varinfo__, -Inf) + Turing.@addlogprob! (; loglikelihood = -Inf) return nothing end for i in eachindex(datakeys) diff --git a/src/ensemble.jl b/src/ensemble.jl index a4480d0..77e6e99 100644 --- a/src/ensemble.jl +++ b/src/ensemble.jl @@ -31,6 +31,23 @@ function ensemble_weights(sol::EnsembleSolution, data_ensem) return weights = predictions \ data end +""" + EnsembleProbForwarder(all_probs) + +Callable used as the `prob_func` of the `EnsembleProblem` returned by +[`bayesian_ensemble`](@ref). It selects the per-trajectory problem from the stored +`all_probs` vector. It supports both the `prob_func(prob, ctx)` interface of newer +SciMLBase (selecting via `ctx.sim_id`) and the legacy `prob_func(prob, i, repeat)` +interface (selecting via the integer index). Storing `all_probs` lets callers recover +the number of trajectories via `enprob.prob_func.all_probs`. +""" +struct EnsembleProbForwarder{P} + all_probs::P +end + +(f::EnsembleProbForwarder)(prob, i::Integer, repeat) = f.all_probs[i] +(f::EnsembleProbForwarder)(prob, ctx) = f.all_probs[ctx.sim_id] + function bayesian_ensemble( probs, ps, datas; noise_prior = InverseGamma(2, 3), @@ -57,6 +74,6 @@ function bayesian_ensemble( @info "$(length(all_probs)) total models" return enprob = EnsembleProblem( - all_probs[1]; prob_func = (prob, ctx) -> all_probs[ctx.sim_id] + all_probs[1]; prob_func = EnsembleProbForwarder(all_probs) ) end diff --git a/test/ensemble.jl b/test/ensemble.jl index 39f9da3..7f1749b 100644 --- a/test/ensemble.jl +++ b/test/ensemble.jl @@ -53,7 +53,9 @@ eqs = [ @mtkbuild sys3 = ODESystem(eqs, t) prob3 = ODEProblem(sys3, [], tspan); probs = [prob, prob2, prob3] -enprob = EnsembleProblem(probs[1]; prob_func = (prob, ctx) -> probs[ctx.sim_id]) +prob_func(prob, i::Integer, repeat) = probs[i] +prob_func(prob, ctx) = probs[ctx.sim_id] +enprob = EnsembleProblem(probs[1]; prob_func = prob_func) sol = solve(enprob, Tsit5(); saveat = 1, trajectories = length(probs)); From b419b8a445bc8c3649bd0803c1a469d1f21ac211 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 20 Jun 2026 09:21:07 -0400 Subject: [PATCH 3/5] Make bayesian_datafit chain extraction + sample rejection cross-version (Downgrade) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous commit fixed the Core/Datafit reds on the current stack but two bayesian_datafit details still broke the Downgrade group (Turing 0.35 / MCMCChains 6 / DynamicPPL 0.30): 1. Chain sample extraction. `chain[@varname(pprior[i])]` works on newer Turing's FlexiChains.VNChain but errors on the legacy MCMCChains.Chains: MethodError: no method matching getindex(::MCMCChains.Chains{...}, ::VarName{:pprior, IndexLens{Tuple{Int64}}}) Added `_pprior_samples(chain, i)` which tries the VarName form and falls back to the legacy `chain["pprior[i]"]` string key, supporting both backends. 2. Sample rejection. `@addlogprob! (; loglikelihood = -Inf)` (NamedTuple form) only exists on DynamicPPL 0.41+. Switched to the scalar `@addlogprob! -Inf`, which is valid on both 0.30 (added to the log-prob) and 0.41 (routed to the LogLikelihood accumulator), and still avoids the removed `acclogp!!(__varinfo__, ::Float64)`. Verified the VarName→string fallback against a real MCMCChains.Chains (the exact type the Downgrade run reported), and re-verified bayesian_datafit on the current stack (Julia 1.12, Turing 0.45/FlexiChains): both forms pass the variance-reduction assertion, including runs that actually hit the rejection branch (solves aborting under ForwardDiff), with no acclogp!! error. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- src/datafit.jl | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/datafit.jl b/src/datafit.jl index 8b7208c..12bf499 100644 --- a/src/datafit.jl +++ b/src/datafit.jl @@ -1,3 +1,23 @@ +""" + _pprior_samples(chain, i) + +Extract the posterior samples of `pprior[i]` from a Turing sampling result, supporting +both chain backends across Turing versions: the `FlexiChains.VNChain` returned by newer +Turing (indexed by `@varname(pprior[i])`) and the legacy `MCMCChains.Chains` +(indexed by the `"pprior[i]"` string key). +""" +function _pprior_samples(chain, i) + vn = @varname(pprior[i]) + samples = try + chain[vn] + catch err + err isa Union{MethodError, KeyError, ArgumentError} || + rethrow(err) + chain["pprior[" * string(i) * "]"] + end + return collect(samples)[:] +end + function l2loss(pvals, (prob, pkeys, t, data)::Tuple{Vararg{Any, 4}}) p = Pair.(pkeys, pvals) prob = remake(prob, tspan = (prob.tspan[1], t[end]), p = p) @@ -223,7 +243,7 @@ Turing.@model function bayesianODE(prob, t, pdist, pkeys, data, noise_prior) prob = remake(prob, tspan = (prob.tspan[1], t[end]), p = Pair.(pkeys, pprior)) sol = solve(prob, saveat = t) if !SciMLBase.successful_retcode(sol) - Turing.@addlogprob! (; loglikelihood = -Inf) + Turing.@addlogprob! -Inf return nothing end for i in eachindex(data) @@ -249,7 +269,7 @@ Turing.@model function bayesianODE( prob = remake(prob, tspan = (prob.tspan[1], lastt), p = Pair.(pkeys, pprior)) sol = solve(prob) if !SciMLBase.successful_retcode(sol) - Turing.@addlogprob! (; loglikelihood = -Inf) + Turing.@addlogprob! -Inf return nothing end for i in eachindex(datakeys) @@ -308,7 +328,7 @@ function bayesian_datafit( progress = true ) return [ - Pair(p[i].first, collect(chain[@varname(pprior[i])])[:]) + Pair(p[i].first, _pprior_samples(chain, i)) for i in eachindex(p) ] end @@ -333,7 +353,7 @@ function bayesian_datafit( progress = true ) return [ - Pair(p[i].first, collect(chain[@varname(pprior[i])])[:]) + Pair(p[i].first, _pprior_samples(chain, i)) for i in eachindex(p) ] end From e3d26bba83be273ffcbdbad3e85701640e87bb67 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 20 Jun 2026 10:53:58 -0400 Subject: [PATCH 4/5] Fix get_sensitivity EnsembleProblem prob_func arity (Core sensitivity.jl) `_get_sensitivity`'s internal `prob_func(prob, i, repeat)` only defined the legacy 3-arg form, so on current SciMLBase (3.x) the EnsembleProblem solve called it with the 2-arg `(prob, ctx)` form and errored: MethodError: no method matching (::#prob_func#25)(::ODEProblem, ::EnsembleContext) Added the `(prob, ctx) -> remake(...; p = ...[:, ctx.sim_id])` method alongside the integer-index one (same cross-version pattern as the ensemble fix), and read per-trajectory solutions via `sol.u[i]` (matching the EnsembleSolution access used in src/ensemble.jl). This is the same upstream EnsembleProblem prob_func API change; it surfaced in the Core group only after the ensemble.jl fix let Core run past it. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- src/sensitivity.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sensitivity.jl b/src/sensitivity.jl index ea52199..4f5f0ed 100644 --- a/src/sensitivity.jl +++ b/src/sensitivity.jl @@ -2,7 +2,8 @@ function _get_sensitivity(prob, t, x, pbounds; samples) boundvals = getfield.(pbounds, :second) boundkeys = getfield.(pbounds, :first) f = function (p) - prob_func(prob, i, repeat) = remake(prob; p = Pair.(boundkeys, p[:, i])) + prob_func(prob, i::Integer, repeat) = remake(prob; p = Pair.(boundkeys, p[:, i])) + prob_func(prob, ctx) = remake(prob; p = Pair.(boundkeys, p[:, ctx.sim_id])) ensemble_prob = EnsembleProblem(prob, prob_func = prob_func) sol = solve( ensemble_prob, nothing, EnsembleThreads(); saveat = t, @@ -11,11 +12,11 @@ function _get_sensitivity(prob, t, x, pbounds; samples) out = zeros(size(p, 2)) if x isa Function for i in 1:size(p, 2) - out[i] = x(sol[i]) + out[i] = x(sol.u[i]) end else for i in 1:size(p, 2) - out[i] = sol[i](t; idxs = x) + out[i] = sol.u[i](t; idxs = x) end end return out From de5f31095edd8743ed343c55b1385305d9c875ad Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 22 Jun 2026 23:46:47 -0400 Subject: [PATCH 5/5] Fix prob_violating_threshold for Symbolics 7 (Core threshold.jl errors) Continuing this branch's "adapt to upstream SciML churn" theme, the Core group's threshold.jl surfaced two more errors once the ensemble/sensitivity fixes let Core run past them. `prob_violating_threshold`'s inner `g(sol, p)` introspected each symbolic threshold inequality via `threshold.val.f` and `threshold.val.arguments`. On Symbolics 7 the `BasicSymbolic` representation is backed by Moshi `@data`, so direct `.f`/`.arguments` field access errors with `unknown field name: arguments` (from Moshi's getproperty), failing both `prob_violating_threshold` tests (threshold.jl:158 `> 0.99` and :165 `< 0.01`) with that error. Replaced the field access with the public `operation`/`arguments` accessors via a small `_threshold_violation` helper. Symbolics also canonicalizes comparisons (`x > 10.0` is stored as `<(10.0, x)`), moving the constant to either operand, so the helper detects which argument is the numeric bound and reconstructs the violating side (upper: `maximum(state) > bound`; lower: `minimum(state) < bound`) accordingly. Verified locally on Julia 1.12 (SciMLBase 3.22, ModelingToolkit 11.28, Symbolics 7): test/threshold.jl now reports 10 passed / 1 failed / 0 errored (was 8/1/2). The two `unknown field name: arguments` errors are gone; both `prob_violating_threshold` assertions pass. The remaining single failure (threshold.jl:123, `2.97 < 2`) is a pre-existing, machine-load-dependent flake in `optimal_parameter_intervention_for_threshold`, which optimizes with NLopt's stochastic `GN_ISRES` under a wall-clock `maxtime=60` budget and no fixed RNG seed; it is unrelated to these changes and is not touched here. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- src/threshold.jl | 35 +++++++++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/src/threshold.jl b/src/threshold.jl index 3ab0039..9a22f0e 100644 --- a/src/threshold.jl +++ b/src/threshold.jl @@ -32,6 +32,30 @@ function get_threshold(prob, obs, threshold; alg = nothing, kw...) return sol.t[end] end +# Decompose a symbolic threshold inequality (e.g. `x > 10.0`) into the state it +# constrains, the numeric bound, and whether the violating side is the upper one +# (`maximum(state) > bound`) or the lower one (`minimum(state) < bound`). +# Symbolics canonicalizes comparisons, so `x > 10.0` is stored as `<(10.0, x)`: +# the constant can land on either side of the operator, which this normalizes. +function _threshold_violation(threshold) + v = ModelingToolkit.value(threshold) + op = operation(v) + args = arguments(v) + isconst(z) = ModelingToolkit.value(z) isa Number + if isconst(args[1]) + bound = ModelingToolkit.value(args[1]) + state = args[2] + # `bound op state`: `bound < state` ⟺ `state > bound` (upper violation). + upper = (op === <) || (op === <=) + else + bound = ModelingToolkit.value(args[2]) + state = args[1] + # `state op bound`: `state > bound`/`state >= bound` is the upper violation. + upper = (op === >) || (op === >=) + end + return state, bound, upper +end + """ prob_violating_thresholdd(prob, p, thresholds) @@ -46,16 +70,15 @@ function prob_violating_threshold(prob, p, thresholds) h(x, u, p) = u, remake(prob, p = Pair.(pkeys, [x...])).p # remake does not work well with static arrays function g(sol, p) for threshold in thresholds - if (threshold.val.f == >) || (threshold.val.f == >=) - if maximum(sol[threshold.val.arguments[1]]) > threshold.val.arguments[2] + state, bound, upper = _threshold_violation(threshold) + if upper + if maximum(sol[state]) > bound return 1.0 end - elseif (threshold.val.f == <) || (threshold.val.f == <=) - if minimum(sol[threshold.val.arguments[1]]) < threshold.val.arguments[2] + else + if minimum(sol[state]) < bound return 1.0 end - else - error() end end return 0.0