Add params support to Contractor and Separator

.gitignore

@@ -2,4 +2,5 @@
 *.jl.*.cov
 *.jl.mem
 docs/build/
-Manifest.toml
+Manifest.toml
+papers/

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.

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"

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")

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