diff --git a/Project.toml b/Project.toml index 9f83f9a..37ad6e9 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DiscreteMeasures = "7766d772-2108-41ee-a4bd-11c51440a39b" OUQBase = "01930cae-99d2-7439-8f4f-ace2ece9f1b9" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [sources] CanonicalMoments = { path = "lib/CanonicalMoments" } @@ -24,6 +25,7 @@ Pkg = "1.10" PrecompileTools = "1.2" Reexport = "1.2.2" SafeTestsets = "0.1" +Symbolics = "7.37" Test = "1.10" SciMLTesting = "2" diff --git a/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl b/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl index 62ae7ca..f42ef52 100644 --- a/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl +++ b/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl @@ -2,10 +2,9 @@ module IntervalArithmeticExt using DiscreteMeasures, IntervalArithmetic -println("Load Ext") function DiscreteMeasures.clamp_domain(x::Interval, lb, ub) bounds = interval(lb, ub) - return intersect(x, bounds) + return intersect_interval(x, bounds) end function DiscreteMeasures.clamp_weight(w::Interval{T}, maxw::T = one(T)) where {T} diff --git a/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl b/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl index 14a48ab..ac42b98 100644 --- a/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl +++ b/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl @@ -91,11 +91,12 @@ end end @testset "Interval Ext" begin - @test clamp_domain(interval(-2, -1), 0, 2) == emptyinterval() - @test clamp_domain(interval(-2, 1), 0, 2) == interval(0, 1) - @test clamp_domain(interval(0.1, 0.5), 0, 2) == interval(0.1, 0.5) - @test clamp_domain(interval(0.5, 5), 0, 2) == interval(0.5, 2) - @test clamp_domain(interval(3, 4), 0, 2) == emptyinterval() + # IntervalArithmetic 1.x: empty intervals compare false under `==` even when equal. + @test isequal_interval(clamp_domain(interval(-2, -1), 0, 2), emptyinterval()) + @test isequal_interval(clamp_domain(interval(-2, 1), 0, 2), interval(0, 1)) + @test isequal_interval(clamp_domain(interval(0.1, 0.5), 0, 2), interval(0.1, 0.5)) + @test isequal_interval(clamp_domain(interval(0.5, 5), 0, 2), interval(0.5, 2)) + @test isequal_interval(clamp_domain(interval(3, 4), 0, 2), emptyinterval()) end end diff --git a/lib/OUQBase/src/OUQBase.jl b/lib/OUQBase/src/OUQBase.jl index 89bc13c..b1f8bd9 100644 --- a/lib/OUQBase/src/OUQBase.jl +++ b/lib/OUQBase/src/OUQBase.jl @@ -1,7 +1,7 @@ module OUQBase using ModelingToolkit: ModelingToolkit, @named, OptimizationSystem, get_variables, getbounds, parameters, structural_simplify, unknowns -using Symbolics: Symbolics, Equation, Inequality, Num, wrap +using Symbolics: Symbolics, Equation, Inequality, Num, wrap, ≲, ≳ using SymbolicUtils: SymbolicUtils, BasicSymbolic, @rule, substitute using SymbolicUtils.Rewriters: Chain using TermInterface: arguments @@ -16,6 +16,9 @@ using JuMP: JuMP, @constraint, @objective, @variable, set_lower_bound, set_name, using SciMLBase: SciMLBase, OptimizationFunction, OptimizationProblem using ADTypes: AbstractADType +const _LEQ_RELATIONAL_OPERATOR = (0 ≲ 0).relational_op +const _GEQ_RELATIONAL_OPERATOR = (0 ≳ 0).relational_op + using CanonicalMoments: CanonicalMoments, DiscreteMeasureTransform1 import CanonicalMoments: RawMomentSequence diff --git a/lib/OUQBase/src/interface.jl b/lib/OUQBase/src/interface.jl index 831dd9c..612ee31 100644 --- a/lib/OUQBase/src/interface.jl +++ b/lib/OUQBase/src/interface.jl @@ -203,18 +203,30 @@ function get_group_name(var, ouq_sys::OUQSystem) return get_group_name(var, ouq_sys.admissible_set) end +function _is_random_variable(ex, admissible_set::AdmissibleSet) + wrapped_ex = wrap(ex) + return any(values(admissible_set.random_variable_map)) do variables + if variables isa AbstractVector + return any(variable -> isequal(wrapped_ex, variable), variables) + end + return isequal(wrapped_ex, variables) + end +end + function get_ordered_group_names( expression, admissible_set::AdmissibleSet; ensure_singleton = true, ensure_all = false, ) + vars = get_variables( + expression; + is_atomic = ex -> _is_random_variable(ex, admissible_set), + ) if ensure_singleton - vars = get_variables(expression) @assert length(vars) == 1 "Expression $expression has multiple random variables $vars. Disable `ensure_singleton` if this is intended." return unique([get_group_name(only(vars), admissible_set)]) else - vars = get_variables(expression) @debug "Expression $expression has multiple random variables $vars" unordered_group_names = Set([get_group_name(var, admissible_set) for var in vars]) if ensure_all diff --git a/lib/OUQBase/src/operators.jl b/lib/OUQBase/src/operators.jl index adb9562..21d5659 100644 --- a/lib/OUQBase/src/operators.jl +++ b/lib/OUQBase/src/operators.jl @@ -1,14 +1,20 @@ -import Symbolics: Operator, value -import SymbolicUtils: Term, symtype +import Symbolics: value +import SymbolicUtils: term -struct 𝔼_ <: Operator end +struct 𝔼_ end const 𝔼 = 𝔼_() # To get same object -(::𝔼_)(x) = Term{symtype(x)}(𝔼, Any[x]) +# `term` (rather than constructing `Term{T}` directly) picks the SymVariant +# type parameter itself and infers the symbolic type via `promote_symtype`. +(::𝔼_)(x) = term(𝔼, x) (::𝔼_)(x::Num) = Num(𝔼(value(x))) SymbolicUtils.promote_symtype(::𝔼_, x) = x +# Disambiguates against SymbolicUtils's generic `promote_symtype(::Operator, +# ::Type{T}) where T`, which is equally specific to the untyped method above +# when called with a `Type` argument (e.g. `symtype(x) === Real`). +SymbolicUtils.promote_symtype(::𝔼_, ::Type{T}) where {T} = T Base.nameof(::𝔼_) = :𝔼 SymbolicUtils.isbinop(::𝔼_) = false @@ -24,14 +30,15 @@ struct ℙ <: Operator end =# -struct ℙ_ <: Operator end +struct ℙ_ end const ℙ = ℙ_() # To get same object -(::ℙ_)(x) = Term{symtype(x)}(ℙ, Any[x]) +(::ℙ_)(x) = term(ℙ, x) (::ℙ_)(x::Num) = Num(ℙ(value(x))) SymbolicUtils.promote_symtype(::ℙ_, x) = x +SymbolicUtils.promote_symtype(::ℙ_, ::Type{T}) where {T} = T Base.nameof(::ℙ_) = :ℙ SymbolicUtils.isbinop(::ℙ_) = false diff --git a/lib/OUQBase/src/reduction_transformations/canonical_moments.jl b/lib/OUQBase/src/reduction_transformations/canonical_moments.jl index 1ea4761..c9571a7 100644 --- a/lib/OUQBase/src/reduction_transformations/canonical_moments.jl +++ b/lib/OUQBase/src/reduction_transformations/canonical_moments.jl @@ -6,7 +6,7 @@ function get_raw_moment_order(equation::Union{Equation, Inequality}, random_var: if Symbolics.wrap(_eq.lhs) === random_var return 1 end - order = _eq.lhs + order = Symbolics.value(_eq.lhs) if !isa(order, Int64) error( "Equation $(equation) is not a raw moment equation of the form: 𝔼(Q^n) ~ where n is an Integer", @@ -21,8 +21,9 @@ function build_raw_moment_sequence( ) num_group_cons = length(constraints) raw_moment_sequence = fill(NaN, num_group_cons) - for (i, constraint) in enumerate(constraints) - raw_moment_sequence[get_raw_moment_order(constraint, random_var)] = constraint.rhs + for constraint in constraints + moment_order = get_raw_moment_order(constraint, random_var) + raw_moment_sequence[moment_order] = Symbolics.value(constraint.rhs) end !any(isnan, raw_moment_sequence) || error("Raw moments have holes") lb, ub = getbounds(random_var) @@ -41,6 +42,17 @@ function create_raw_moments_map(admissible_set::AbstractAdmissibleSet) return raw_moments_map end +function _evaluate_condition(condition::Inequality, substitutions) + lhs = Symbolics.value(substitute(condition.lhs, substitutions; fold = Val(true))) + rhs = Symbolics.value(substitute(condition.rhs, substitutions; fold = Val(true))) + if condition.relational_op == _LEQ_RELATIONAL_OPERATOR + return lhs <= rhs + elseif condition.relational_op == _GEQ_RELATIONAL_OPERATOR + return lhs >= rhs + end + throw(ArgumentError("Unsupported inequality: $condition")) +end + # This creates the optimization decision variables for the canonical moment problem. # Some code repetition with create_discrete_measure but keeping both makes sense because of the need to provide correct bounds for supports in discrete_measure case. function create_free_variable_vec(group_name::Symbol, n_supports::Int64) @@ -250,7 +262,7 @@ function construct_optimization_problem( if isa(ouq_sys.objective, ProbabilityObjective) # TODO: Dispatch on Probability Function approximation here. ouq_obj_f = - (rand_var_vec) -> Symbolics.evaluate( + (rand_var_vec) -> _evaluate_condition( condition, Dict(constituent_random_variables .=> rand_var_vec), ) diff --git a/lib/OUQBase/src/reduction_transformations/winkler_extremal_measures.jl b/lib/OUQBase/src/reduction_transformations/winkler_extremal_measures.jl index 8e499bf..96b8f7e 100644 --- a/lib/OUQBase/src/reduction_transformations/winkler_extremal_measures.jl +++ b/lib/OUQBase/src/reduction_transformations/winkler_extremal_measures.jl @@ -5,9 +5,9 @@ function convert_inequality_to_jump_leq_lhs( complement = false, tol = 1.0e-10, ) - if ineq.relational_op == Symbolics.leq + if ineq.relational_op == _LEQ_RELATIONAL_OPERATOR jump_leq_lhs = complement ? ineq.rhs - ineq.lhs + tol : ineq.lhs - ineq.rhs - elseif ineq.relational_op == Symbolics.geq + elseif ineq.relational_op == _GEQ_RELATIONAL_OPERATOR jump_leq_lhs = complement ? ineq.lhs - ineq.rhs - tol : ineq.rhs - ineq.lhs else error("Unsupported inequality: $ineq") @@ -386,10 +386,9 @@ function construct_optimization_problem( ensure_all = true, ) - ## In the exact case, we are using the Symbolic.evaluate function # TODO: Dispatch on Probability Function approximation here. reduced_objective = expectation( - (single_support) -> Symbolics.evaluate( + (single_support) -> _evaluate_condition( condition, Dict(constituent_random_variables .=> single_support), ), diff --git a/lib/OUQBase/test/Core/FloodProblem/Q_only.jl b/lib/OUQBase/test/Core/FloodProblem/Q_only.jl index 1a47e5e..4438eb0 100644 --- a/lib/OUQBase/test/Core/FloodProblem/Q_only.jl +++ b/lib/OUQBase/test/Core/FloodProblem/Q_only.jl @@ -1,5 +1,6 @@ using Symbolics, ModelingToolkit using OUQBase +using CanonicalMoments: moments using Test using OptimizationBBO @@ -7,6 +8,10 @@ using OptimizationBBO rand_vars = @random_variables begin Independent(Q, bounds = (160.0, 3580.0)) end +@test OUQBase._evaluate_condition(Q ≲ 160.0, Dict(Q => 160.0)) +@test OUQBase._evaluate_condition(Q ≳ 160.0, Dict(Q => 160.0)) +@test !OUQBase._evaluate_condition(Q ≲ 160.0, Dict(Q => 161.0)) +@test !OUQBase._evaluate_condition(Q ≳ 160.0, Dict(Q => 159.0)) constraints = [𝔼(Q) ~ 1320.42] @@ -41,6 +46,17 @@ ouq_sys_expectation_canonical_moments = OUQSystem(; reduction_alg = StengerCanonicalMoments(), parameters = pars, ) +@test Symbolics.value.( + moments( + only( + values( + raw_moments_map( + ouq_sys_expectation_canonical_moments, + ) + ) + ) + ) +) == [1320.42] ouq_sys_expectation_canonical_moments_analytic = OUQSystem(; objective = objective_expectation, admissible_set, diff --git a/lib/OUQBase/test/qa/qa.jl b/lib/OUQBase/test/qa/qa.jl index b9d53e3..e508c23 100644 --- a/lib/OUQBase/test/qa/qa.jl +++ b/lib/OUQBase/test/qa/qa.jl @@ -18,13 +18,12 @@ run_qa( # Non-public qualified accesses into upstream packages; still non-public in the # resolved upstream majors: # :BasicSymbolic/:isbinop/:promote_symtype - SymbolicUtils - # :evaluate/:geq/:leq - Symbolics # :getdefault - ModelingToolkit # :NoAD/:NullParameters - SciMLBase all_qualified_accesses_are_public = (; ignore = ( - :BasicSymbolic, :NoAD, :NullParameters, :evaluate, :geq, :getdefault, - :isbinop, :leq, :promote_symtype, + :BasicSymbolic, :NoAD, :NullParameters, :getdefault, :isbinop, + :promote_symtype, ), ), ), diff --git a/src/OptimalUncertaintyQuantification.jl b/src/OptimalUncertaintyQuantification.jl index bc5b3ad..58a8d61 100644 --- a/src/OptimalUncertaintyQuantification.jl +++ b/src/OptimalUncertaintyQuantification.jl @@ -3,6 +3,11 @@ module OptimalUncertaintyQuantification using Reexport: @reexport @reexport using OUQBase +# `@random_variables` expands to code that references `Symbolics.@variables` +# in the caller's scope (via `esc`), so `Symbolics` must be bound here for the +# precompile workload below to expand successfully. +using Symbolics: Symbolics + using PrecompileTools: @compile_workload, @setup_workload @setup_workload begin @@ -11,7 +16,7 @@ using PrecompileTools: @compile_workload, @setup_workload Independent(Q, bounds = (0.0, 1.0)) end admissible_set = AdmissibleSet(random_vars, [𝔼(Q) ~ 0.5]) - OUQSystem( + OUQSystem(; objective = 𝔼(Q), admissible_set, reduction_alg = WinklerExtremalMeasures(),