diff --git a/.gitignore b/.gitignore index 44a110c..fa8b6af 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ *.jl.*.cov *.jl.mem docs/build/ -Manifest.toml \ No newline at end of file +Manifest.toml +papers/ diff --git a/CLAUDE.md b/CLAUDE.md index 395e1fc..062f074 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -43,3 +43,66 @@ This file tracks decisions and context discovered while working on this package. - Chained comparisons like `0 <= x^2+y^2 <= 1` don't work (Julia lowers to `&&` which requires `Bool`). Users should use `x^2+y^2 ∈ interval(0, 1)` instead. A `@constraint` macro could fix this — see `future.md`. - ReversePropagation emits many "Method definition overwritten" warnings with Symbolics v7. Harmless but noisy — see `future.md`. + +## ReversePropagation SSAFunction / RuntimeGeneratedFunctions migration (2026-04-20) + +### Versions + +Bumped `ReversePropagation` compat from `"0.3 - 0.4"` to `"0.8"`. Added a direct dependency on `RuntimeGeneratedFunctions = "0.5"` (transitively available through RP, but worth owning directly since we now call `@RuntimeGeneratedFunction` ourselves). Bumped the package version to `0.16.0`. + +### Why + +RP 0.5+ reshapes `cse_equations`, `forward_backward_code`, `forward_backward_contractor`, and `gradient_code` around a new `SSAFunction` value: SSA `code`, an `output` variable, and a `variables` NamedTuple carrying extras like `constraint`. It also replaces `eval` with `@RuntimeGeneratedFunction` so a generated function can be built and called within the same dynamic extent (no world-age error) and with a 4–15× faster build time. + +### Code changes + +**`src/IntervalConstraintProgramming.jl`** — bring in `SSAFunction` and `cse_equations` from RP; call `RuntimeGeneratedFunctions.init(@__MODULE__)` so `@RuntimeGeneratedFunction` works inside this module; export `SSAFunction`. + +**`src/utils.jl`** — `make_function` now wraps `build_function(...)` in `@RuntimeGeneratedFunction` instead of `eval`. Same reasons as upstream RP: no world-age hazard, faster build. + +**`src/contractor.jl`** — `Contractor` now stores the `SSAFunction` it was built from: +- `Contractor(ex, vars)` calls `cse_equations(ex)` once and forwards to +- `Contractor(ssa::SSAFunction, ex, vars)` which calls `forward_backward_contractor(ssa, vars)` (the SSA-taking method). +- `Separator(ex, vars, constraint::Interval)` runs `cse_equations(ex)` once and shares the SSA between `make_function(ex, vars)` (forward evaluator) and `Contractor(ssa, ex, vars)` (HC4Revise), so each constructor does CSE on `ex` only once rather than twice. + +### Known issues + +- This change requires ReversePropagation `0.8.0`, which is not registered yet at the time of writing — to develop/test against it, `]dev` the local RP checkout. + +## Suppressing the "NG" flag via `exact=true` (2026-04-20) + +### Problem + +`IntervalArithmetic` v1 sets a "not guaranteed" (NG) flag on any interval that has ever touched a plain (non-`ExactReal`) `Real`. ICP's generated code is full of those — both the literals from the user's expression (the `1` in `x^2 + y^2 == 1`) and the integer exponents that CSE inserts (the `2` in `x^2`). So even `C.f([3..4, 5..6])` came back as `[33.0, 51.0]_com_NG`. + +### Why not on the symbolic side + +Wrapping literals in the user-facing symbolic expression (`x^2 + y^2 - exact(1)`) does not survive — `Symbolics` normalises `ExactReal` back to a plain `Int`/`Float` during its own simplification. Any fix has to happen on the generated `Expr`, after `Symbolics.build_function` / `ReversePropagation.forward_backward_expr`. + +### Design + +Added an `exactify(ex)` helper in `src/utils.jl` that walks an `Expr` and wraps every non-`Bool` `Real` literal with `IntervalArithmetic.exact(...)`. Threaded an `exact::Bool = false` keyword argument through `Contractor`, `Separator`, and `separator`/`constraint`: + +- `make_function(ex, vars; exact=false)` — runs `exactify` on the `build_function` output when `exact=true`, before `@RuntimeGeneratedFunction`. +- `Contractor(ssa, ex, vars; exact=false)` — when `exact=true`, delegates to `build_fb_contractor`, which replicates RP's outer wrapper locally so it can `exactify` and run `preserve_decorations` on the body before RGF. When `exact=false`, still forwards to `forward_backward_contractor` unchanged (no overhead in the default path). +- `Separator(...; exact=false)` and `separator(...; exact=false)` — threads the keyword through. + +### Decoration preservation through contraction + +IEEE 1788 weakens `intersect_interval(a, b)` to the `trv` decoration because, in general, the result may not be a subset of the range of any function. Every reverse op in `IntervalContractors` (`plus_rev`, `power_rev`, …) is built on top of that same intersect, so both the outer `intersect_interval(_value, _constraint)` that RP emits *and* every narrowing step in the reverse sweep drag `trv` through the result. HC4Revise however is only ever narrowing *within* an already-evaluated function range, so the right answer is `min(decoration(inputs)...)` — not `trv`. + +`preserve_decorations(ex)` walks the body `Expr` returned by `ReversePropagation.forward_backward_expr` and wraps each `IntervalArithmetic.intersect_interval` call and each `ReversePropagation._rev_*` reverse-op call so its output interval(s) are re-decorated with `min(decoration(inputs)...)`. Empty intervals stay `trv` (empty semantics). Because Symbolics' `toexpr` emits the call head as the actual function *value* (not `Expr(:., Module, :name)`), the matchers compare `ex.args[1] === IntervalArithmetic.intersect_interval` and `parentmodule(f) === ReversePropagation && startswith(nameof(f), "_rev_")` directly. + +Re-decoration strategy: `_redecorate(x, d) = interval(inf(x), sup(x), d)` unconditionally replaces the decoration; it does *not* `min` with `decoration(x)`, because the op's own output decoration is exactly the weakened `trv` we're trying to override. Narrowing is a subset, so the pre-op input decoration is the valid upper bound regardless of what the op stamped on. Empty intervals keep their existing (`trv`) decoration. + +### Usage + +```julia +julia> C = constraint(x^2 + y^2 == 1, vars; exact=true); +julia> C.f([3..4, 5..6]) +[33.0, 51.0]_com # no NG +julia> C.contractor((-10..10, -10..10)) +[-1.0, 1.0]_com² # no NG, no trv +``` + +Default behaviour (`exact=false`) is unchanged. Note: the outer of the full `Separator` triple can still end up `_trv` because it is computed via `boundary ⊔ corner` in `(SS::Separator)(X)` — hull (`⊔`) intrinsically weakens in IEEE 1788, and that is the intended semantics. diff --git a/Project.toml b/Project.toml index 56e51d3..f547d81 100644 --- a/Project.toml +++ b/Project.toml @@ -1,20 +1,20 @@ name = "IntervalConstraintProgramming" uuid = "138f1668-1576-5ad7-91b9-7425abbf3153" -version = "0.15.0" +version = "0.16.0" [deps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" IntervalBoxes = "43d83c95-ebbb-40ec-8188-24586a1458ed" -IntervalContractors = "15111844-de3b-5229-b4ba-526f2f385dc9" ReversePropagation = "527681c1-8309-4d3f-8790-caf822a419ba" +RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] IntervalArithmetic = "1" IntervalBoxes = "0.3" -IntervalContractors = "0.6" -ReversePropagation = "0.4.1" +ReversePropagation = "0.8" +RuntimeGeneratedFunctions = "0.5" StaticArrays = "1" Symbolics = "6 - 7" julia = "1.10" diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl index 12b02be..b727196 100644 --- a/src/IntervalConstraintProgramming.jl +++ b/src/IntervalConstraintProgramming.jl @@ -3,7 +3,9 @@ __precompile__() module IntervalConstraintProgramming using IntervalArithmetic, IntervalArithmetic.Symbols -using IntervalContractors +## IntervalContractors removed on this branch — see ReversePropagation +## branch lorentz-no-icontractors. Was previously imported only to bind +## `reverse_operations`, which was never referenced inside this package. using IntervalBoxes using Symbolics @@ -14,6 +16,10 @@ using Symbolics: @register_symbolic using StaticArrays using ReversePropagation +using ReversePropagation: SSAFunction, cse_equations + +using RuntimeGeneratedFunctions +RuntimeGeneratedFunctions.init(@__MODULE__) # using MacroTools @@ -50,11 +56,11 @@ export add_constraint!, variables, add_variables!, - ConstraintProblem + ConstraintProblem, + SSAFunction export ¬ -const reverse_operations = IntervalContractors.reverse_operations include("utils.jl") diff --git a/src/contractor.jl b/src/contractor.jl index 3181bbd..891a4d8 100644 --- a/src/contractor.jl +++ b/src/contractor.jl @@ -1,26 +1,110 @@ abstract type AbstractContractor end -struct Contractor{V, E, CC} <: AbstractContractor +""" + Contractor(ex, vars; exact=false, params=[]) + Contractor(ssa::SSAFunction, ex, vars; exact=false, params=[]) + +HC4Revise forward--backward contractor for the constraint `ex(vars...; params...) ∈ constraint`. + +The symbolic expression is first lowered to SSA form via +`ReversePropagation.cse_equations`; the resulting `SSAFunction` is stored on the +contractor so it can be reused (e.g. shared between a `Separator` and other +contractors over the same expression). + +If `exact=true`, every numeric literal in the generated forward--backward code +is wrapped in `IntervalArithmetic.exact` (see [`exactify`](@ref)) and every +`intersect_interval` / `_rev_*` call site is wrapped so its output decoration +is `min(decoration(inputs)...)` instead of the IEEE 1788 default `trv` (see +[`preserve_decorations`](@ref)). Together these suppress the "NG" flag and +preserve `com` / `dac` decorations through contraction. + +If `params` is non-empty, the contractor compiles a function of the form +`(__inputs, __constraint, __params) -> ...`. At call time, supply `param_vals` +as a positional third argument: `CC(X, constraint, param_vals)`. This lets the +same compiled contractor be re-used at runtime with different parameter values +— useful when a separator must be applied many times with different constants +(e.g. once per blocker disc), avoiding rebuild costs that scale with the number +of distinct parameter values. +""" +struct Contractor{V, E, S, CC, P} <: AbstractContractor vars::V ex::E + ssa::S contractor::CC + params::P end -Contractor(ex, vars) = Contractor(vars, ex, forward_backward_contractor(ex, vars)) +Contractor(ex, vars; exact::Bool = false, params = []) = + Contractor(cse_equations(ex), ex, vars; exact = exact, params = params) + +Contractor(ssa::SSAFunction, ex, vars; exact::Bool = false, params = []) = + Contractor(vars, ex, ssa, + build_fb_contractor(ssa, vars; exact = exact, params = params), + params) + +# Build the forward--backward contractor function. When `exact=false`, just +# forward to RP's `forward_backward_contractor`. When `exact=true`, re-assemble +# the outer wrapper so we can (a) mark every Real literal as +# `IntervalArithmetic.exact` (so constraint literals and CSE-inserted exponents +# don't drag the NG flag through the result) and (b) clamp decorations through +# intersect/reverse-op call sites to `min(input decorations)` instead of the +# IEEE 1788 default `trv` weakening. +# +# `params` (if non-empty) becomes a third positional argument `__params` to +# the generated function — see RP's `forward_backward_contractor` for the +# wiring of the same convention in the non-exact path. +function build_fb_contractor(ssa::SSAFunction, vars; exact::Bool = false, params = []) + exact || return forward_backward_contractor(ssa, vars, params) + + code, final_var, constraint_var = + ReversePropagation.forward_backward_expr(ssa, vars, params) + + input_vars = toexpr(Symbolics.MakeTuple(vars)) + final = toexpr(final_var) + constraint = toexpr(constraint_var) + + body = preserve_decorations(exactify(code)) + + if isempty(params) + function_code = :( + (__inputs, __constraint) -> begin + $input_vars = __inputs + $constraint = __constraint + $body + return $input_vars, $(final) + end + ) + else + params_tuple = toexpr(Symbolics.MakeTuple(params)) + function_code = :( + (__inputs, __constraint, __params) -> begin + $input_vars = __inputs + $constraint = __constraint + $params_tuple = __params + $body + return $input_vars, $(final) + end + ) + end + + return @RuntimeGeneratedFunction(function_code) +end (CC::Contractor)(X, constraint=interval(0.0)) = IntervalBox(CC.contractor(X, constraint)[1]) +(CC::Contractor)(X, constraint, param_vals::Tuple) = IntervalBox(CC.contractor(X, constraint, param_vals)[1]) abstract type AbstractSeparator end "A separator models the inside and outside of a constraint set using a forward--backward contractor" -struct Separator{V,E,C,F,R} <: AbstractSeparator +struct Separator{V,E,C,F,R,P} <: AbstractSeparator vars::V ex::E constraint::C f::F contractor::R + params::P end # Base.show(io::IO, S::Separator) = print(io, "Separator($(S.ex) ∈ $(S.constraint), vars = $(join(S.vars, ", ")))") @@ -33,35 +117,45 @@ function Base.show(io::IO, S::AbstractSeparator) end end -function Separator(orig_expr, vars) +function Separator(orig_expr, vars; exact::Bool = false, params = []) ex, constraint = analyse(orig_expr) - return Separator(ex, vars, constraint) + return Separator(ex, vars, constraint; exact = exact, params = params) end -Separator(ex, vars, constraint::Interval) = Separator(vars, ex, constraint, make_function(ex, vars), Contractor(ex, vars)) +function Separator(ex, vars, constraint::Interval; exact::Bool = false, params = []) + ssa = cse_equations(ex) + return Separator(vars, ex, constraint, + make_function(ex, vars; exact = exact, params = params), + Contractor(ssa, ex, vars; exact = exact, params = params), + params) +end -function separate_infinite_box(S::Separator, X::IntervalBox) +function separate_infinite_box(S::Separator, X::IntervalBox, param_vals=nothing) # for an box that extends to infinity we cannot evaluate at a corner # so use the old method instead where we do inner and outer contractors separately C = S.contractor a, b = inf(S.constraint), sup(S.constraint) - inner = C(X, interval(a, b)) + contract = param_vals === nothing ? + (cset) -> C(X, cset) : + (cset) -> C(X, cset, param_vals) + + inner = contract(interval(a, b)) # to compute outer, we contract with respect to the complement of `a..b`: local outer if a == -Inf - outer = C(X, interval(b, Inf)) + outer = contract(interval(b, Inf)) elseif b == Inf - outer = C(X, interval(-Inf, a)) + outer = contract(interval(-Inf, a)) else # the complement is a union of two pieces - outer1 = C(X, interval(-Inf, a)) - outer2 = C(X, interval(b, Inf)) + outer1 = contract(interval(-Inf, a)) + outer2 = contract(interval(b, Inf)) outer = outer1 ⊔ outer2 end @@ -73,14 +167,23 @@ end "Returns boundary, inner, outer" -function (SS::Separator)(X) +(SS::Separator)(X) = _separate(SS, X, nothing) + +# Atomic separators that were built with no params silently ignore param_vals, +# so mixed compositions (param + non-param via ⊓/⊔) work without special-casing. +(SS::Separator)(X, param_vals::Tuple) = + _separate(SS, X, isempty(SS.params) ? nothing : param_vals) + +function _separate(SS::Separator, X, param_vals) if any(x -> isinf(diam(x)), X) - return separate_infinite_box(SS, X) + return separate_infinite_box(SS, X, param_vals) end # using the contractor to compute the boundary: - boundary = SS.contractor(X) # contract with respect to 0, which is always the boundary + boundary = param_vals === nothing ? + SS.contractor(X) : # contract w.r.t. 0 + SS.contractor(X, interval(0.0), param_vals) # extend the boundary by evaluating at corners of the box to determine inner and outer: @@ -90,14 +193,18 @@ function (SS::Separator)(X) inner = boundary outer = boundary - lb_image = SS.f(lb) + eval_at = param_vals === nothing ? + (corner) -> SS.f(corner) : + (corner) -> SS.f(corner, param_vals) + + lb_image = eval_at(lb) if !isempty_interval(lb_image) && issubset_interval(lb_image, SS.constraint) inner = inner ⊔ lb else outer = outer ⊔ lb end - ub_image = SS.f(ub) + ub_image = eval_at(ub) if !isempty_interval(ub_image) && issubset_interval(ub_image, SS.constraint) inner = inner ⊔ ub else diff --git a/src/set_operations.jl b/src/set_operations.jl index 4b6e06c..bd565af 100644 --- a/src/set_operations.jl +++ b/src/set_operations.jl @@ -1,8 +1,13 @@ -struct CombinationSeparator{V, E, F} <: AbstractSeparator +# `f` : X -> (boundary, inner, outer) — non-parametric call +# `fp` : (X, params) -> (boundary, inner, outer) — parametric call +# Each combinator builds both, so the same composite separator can be invoked +# either way; atomic operands without params silently ignore `params`. +struct CombinationSeparator{V, E, F, FP} <: AbstractSeparator vars::V ex::E f::F + fp::FP end # TODO: Replace intersection with composition @@ -25,9 +30,23 @@ function ⊓(S1::AbstractSeparator, S2::AbstractSeparator) end + fp = (X, p) -> begin + + boundary1, inner1, outer1 = S1(X, p) + boundary2, inner2, outer2 = S2(X, p) + + inner = inner1 ⊓ inner2 + outer = outer1 ⊔ outer2 + + boundary = inner ⊓ outer + + return (boundary, inner, outer) + + end + ex = (S1.ex) ⊓ (S2.ex) - return CombinationSeparator(vars, ex, f) + return CombinationSeparator(vars, ex, f, fp) end @@ -49,9 +68,23 @@ function ∘(S1::AbstractSeparator, S2::AbstractSeparator) end + fp = (X, p) -> begin + + boundary1, inner1, outer1 = S1(X, p) + boundary2, inner2, outer2 = S2(X, p) + + inner = inner1 ⊓ inner2 + outer = outer1 ⊔ outer2 + + boundary = inner ⊓ outer + + return (boundary, inner, outer) + + end + ex = (S1.ex) ⊓ (S2.ex) - return CombinationSeparator(vars, ex, f) + return CombinationSeparator(vars, ex, f, fp) end @@ -75,9 +108,23 @@ function ⊔(S1::AbstractSeparator, S2::AbstractSeparator) end + fp = (X, p) -> begin + + boundary1, inner1, outer1 = S1(X, p) + boundary2, inner2, outer2 = S2(X, p) + + inner = inner1 ⊔ inner2 + outer = outer1 ⊓ outer2 + + boundary = inner ⊓ outer + + return (boundary, inner, outer) + + end + ex = (S1.ex) ⊔ (S2.ex) - return CombinationSeparator(vars, ex, f) + return CombinationSeparator(vars, ex, f, fp) end ⊓ @@ -96,9 +143,17 @@ function ¬(S::AbstractSeparator) end + fp = (X, p) -> begin + + boundary, inner, outer = S(X, p) + + return (boundary, outer, inner) + + end + ex = ¬(S.ex) - return CombinationSeparator(vars, ex, f) + return CombinationSeparator(vars, ex, f, fp) end @@ -107,6 +162,7 @@ Base.:!(S::AbstractSeparator) = ¬(S) (S::CombinationSeparator)(X) = S.f(X) +(S::CombinationSeparator)(X, p::Tuple) = S.fp(X, p) Base.setdiff(S1::AbstractSeparator, S2::AbstractSeparator) = S1 ⊓ ¬(S2) diff --git a/src/utils.jl b/src/utils.jl index 64b76cd..65a0bd7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,6 +1,123 @@ -make_function(ex, vars) = eval(build_function(ex, vars, nanmath=false)) +""" + exactify(ex) + +Walk `ex` and replace every non-`Bool` `Real` literal with `IntervalArithmetic.exact(literal)`. + +Used on the `Expr` produced by symbolic code generation (`Symbolics.build_function` +or `ReversePropagation.forward_backward_expr`) so that, when the generated function +is called on `Interval` arguments, numeric literals from the user's constraint +(like the `1` in `x^2 + y^2 == 1`) are treated as exact and do not introduce the +"NG" (not-guaranteed) flag in the result. +""" +exactify(ex::Bool) = ex +exactify(ex::Real) = :($(IntervalArithmetic.exact)($ex)) +exactify(ex::Expr) = Expr(ex.head, exactify.(ex.args)...) +exactify(ex) = ex + + +# --- decoration preservation for HC4Revise --------------------------------- +# +# IEEE 1788 weakens `intersect_interval(a, b)` to the `trv` decoration because, +# in general, the intersection may not be the range of any function. In +# HC4Revise however each intersection/reverse op is narrowing *within* an +# already-evaluated function range, so `min(decoration(inputs)...)` is the +# right answer. The helpers below re-decorate the output accordingly; the +# Expr walker `preserve_decorations!` rewrites the call sites that RP's +# `forward_backward_expr` emits so the weakening never takes hold. + +_dec_of(x::Interval) = decoration(x) +_dec_of(x) = IntervalArithmetic.com + +_min_input_decoration(inputs::Tuple) = minimum(_dec_of, inputs) + +# Re-decorate an interval with `d`. Used to undo the `trv` weakening that +# `intersect_interval` and reverse ops perform: in HC4 the output is always a +# subset of one of the inputs, so `d = min(decoration(inputs)...)` is the valid +# decoration, regardless of what the op itself stamped on. Empty intervals +# stay as-is (`trv` / empty semantics). +_redecorate(x::Interval, d::IntervalArithmetic.Decoration) = + isempty_interval(x) ? x : interval(inf(x), sup(x), d) +_redecorate(x, ::IntervalArithmetic.Decoration) = x + +"Wrap an `intersect_interval`-style call: single interval result." +function preserve_dec_intersect(result::Interval, inputs::Tuple) + return _redecorate(result, _min_input_decoration(inputs)) +end + +"Wrap a reverse-op call: tuple result, each interval element re-decorated." +function preserve_dec_reverse(result::Tuple, inputs::Tuple) + d = _min_input_decoration(inputs) + return map(x -> _redecorate(x, d), result) +end + + +# ---- Expr walker ---------------------------------------------------------- +# +# Matches the two forms RP's generated body emits: +# _lhs = (IntervalArithmetic.intersect_interval)(a, b) +# (t1, t2, …) = (ReversePropagation.var"_rev_OP")(arg1, arg2, …) +# and rewrites the rhs so that output decorations are clamped to +# `min(decoration(arg_i))`. + +# Symbolics' `toexpr` resolves call heads to actual function *values* (not +# `Expr(:., Module, QuoteNode(name))`), so these checks compare against the +# function objects directly. + +_is_intersect_call(ex) = false +function _is_intersect_call(ex::Expr) + ex.head === :call || return false + return ex.args[1] === IntervalArithmetic.intersect_interval +end + +_is_reverse_op_call(ex) = false +function _is_reverse_op_call(ex::Expr) + ex.head === :call || return false + f = ex.args[1] + f isa Function || return false + name = string(nameof(f)) + return startswith(name, "_rev_") && + parentmodule(f) === ReversePropagation +end + +""" + preserve_decorations(ex) + +Walk `ex` and wrap calls to `IntervalArithmetic.intersect_interval` and +`ReversePropagation._rev_*` so that each output interval's decoration is +clamped to `min(decoration(inputs)...)` instead of being weakened to `trv`. + +Called on the body `Expr` returned by +`ReversePropagation.forward_backward_expr` before `@RuntimeGeneratedFunction`. +""" +preserve_decorations(ex) = ex +function preserve_decorations(ex::Expr) + if ex.head === :(=) && ex.args[2] isa Expr && _is_intersect_call(ex.args[2]) + lhs = ex.args[1] + call = ex.args[2] + inputs = Expr(:tuple, call.args[2:end]...) + new_rhs = :($(preserve_dec_intersect)($call, $inputs)) + return Expr(:(=), lhs, new_rhs) + elseif ex.head === :(=) && ex.args[2] isa Expr && _is_reverse_op_call(ex.args[2]) + lhs = ex.args[1] + call = ex.args[2] + inputs = Expr(:tuple, call.args[2:end]...) + new_rhs = :($(preserve_dec_reverse)($call, $inputs)) + return Expr(:(=), lhs, new_rhs) + else + return Expr(ex.head, map(preserve_decorations, ex.args)...) + end +end + + +make_function(ex, vars; exact::Bool = false, params = []) = begin + code = isempty(params) ? + build_function(ex, vars, nanmath=false) : + build_function(ex, vars, params, nanmath=false) + exact && (code = exactify(code)) + @RuntimeGeneratedFunction(code) +end """ @@ -41,53 +158,55 @@ end -function separator(ex, vars) - ex2 = ex +function separator(ex, vars; exact::Bool = false, params = []) + ex2 = ex if ex isa Num ex2 = value(ex) end - + op = operation(ex2) if op == ¬ arg = arguments(ex2)[1] - return ¬(separator(arg, vars)) + return ¬(separator(arg, vars; exact = exact, params = params)) end lhs, rhs = arguments(ex2) if op == & - return separator(lhs, vars) ⊓ separator(rhs, vars) + return separator(lhs, vars; exact = exact, params = params) ⊓ + separator(rhs, vars; exact = exact, params = params) elseif op == | - return separator(lhs, vars) ⊔ separator(rhs, vars) + return separator(lhs, vars; exact = exact, params = params) ⊔ + separator(rhs, vars; exact = exact, params = params) elseif op ∈ (≤, <) constraint = interval(-Inf, 0) - Separator(Num(lhs - rhs), vars, constraint) + Separator(Num(lhs - rhs), vars, constraint; exact = exact, params = params) elseif op ∈ (≥, >) constraint = interval(0, +Inf) - Separator(Num(lhs - rhs), vars, constraint) + Separator(Num(lhs - rhs), vars, constraint; exact = exact, params = params) elseif op == (==) constraint = interval(0, 0) - Separator(Num(lhs - rhs), vars, constraint) + Separator(Num(lhs - rhs), vars, constraint; exact = exact, params = params) else - return Separator(ex, vars, interval(0, 0)) # implicit "== 0" + return Separator(ex, vars, interval(0, 0); exact = exact, params = params) # implicit "== 0" end - + end -function separator(ex) +function separator(ex; exact::Bool = false) vars = Symbolics.get_variables(ex) - return separator(ex, vars) + return separator(ex, vars; exact = exact) end diff --git a/test/runtests.jl b/test/runtests.jl index 46d0e32..69cb6c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -116,6 +116,45 @@ end # end +@testset "Parametric Contractor" begin + @variables x, y, r + + # Build a circle-of-radius-r contractor ONCE; evaluate at different r. + C = Contractor(x^2 + y^2, [x, y]; params = [r]) + X = IntervalBox(-Inf..Inf, 2) + + # r = 1: x^2 + y^2 ∈ [-Inf, 1] ⇒ x, y ∈ [-1, 1] + @test eq(C(X, -Inf..1, (interval(1.0),)), IntervalBox(-1..1, 2)) + + # Same compiled contractor at r = 2: x, y ∈ [-2, 2] + @test eq(C(X, -Inf..4, (interval(2.0),)), IntervalBox(-2..2, 2)) +end + +@testset "Parametric Separator" begin + @variables x, y, cx, cy + + II = -100..100 + X = IntervalBox(II, 2) + + # Disc of radius 1 centred at (cx, cy), with (cx, cy) parametric. + S = Separator((x - cx)^2 + (y - cy)^2 <= 1, [x, y]; params = [cx, cy]) + + # Centre (0, 0): inner = [-1, 1]^2. + boundary, inner, outer = S(X, (interval(0.0), interval(0.0))) + @test eq(inner, IntervalBox(-1..1, 2)) + + # Centre (5, 0): inner = [4, 6] × [-1, 1] — same compiled separator. + boundary, inner, outer = S(X, (interval(5.0), interval(0.0))) + @test eq(inner, IntervalBox(4..6, -1..1)) + + # ⊓ with a non-parametric separator: param call should still work, + # the non-parametric leaf ignores the params tuple. + Sx = Separator(x >= 0, [x, y]) # half-plane, no params + Scomb = S ⊓ Sx + boundary, inner, outer = Scomb(X, (interval(0.0), interval(0.0))) + @test eq(inner, IntervalBox(0..1, -1..1)) +end + @testset "Paving a 3D torus" begin vars = @variables x, y, z