From 16aaf50a7514cfb609d7ebcdba82507eeda3a078 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Thu, 2 Apr 2026 03:35:13 -0500 Subject: [PATCH 1/8] Update dependencies to latest versions, fix IntervalArithmetic v1.0 compatibility MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Update compat bounds: IntervalArithmetic 1, IntervalBoxes 0.3, IntervalContractors 0.6, ReversePropagation 0.4, Symbolics 7. IntervalArithmetic v1.0 follows IEEE 1788 and deliberately does not define Base.isequal/Base.hash for Interval. This broke @register_symbolic x ∈ y::Interval since SymbolicUtils needs isequal/hash for hash-consing. Instead of type-pirating those methods, decompose x ∈ interval(a,b) into (x >= a) & (x <= b) at the symbolic level. Also fix pre-existing bug in separator() where & and | used Base.intersect/union instead of ⊓/⊔ (defined for AbstractSeparator in set_operations.jl). Co-Authored-By: Claude Opus 4.6 (1M context) --- CLAUDE.md | 45 ++++++++++++++++++++++++++++ Project.toml | 10 +++---- future.md | 18 +++++++++++ src/IntervalConstraintProgramming.jl | 8 +++-- src/contractor.jl | 10 +++++-- src/utils.jl | 20 ++++++------- 6 files changed, 92 insertions(+), 19 deletions(-) create mode 100644 CLAUDE.md create mode 100644 future.md diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..395e1fc --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,45 @@ +# CLAUDE.md + +This file tracks decisions and context discovered while working on this package. +(This instruction itself: always write discoveries and decisions into CLAUDE.md.) + +## Dependency Update (2026-04-01) + +### Versions updated in Project.toml [compat] + +| Package | Old compat | New compat | Latest version | +|----------------------|-------------|------------|----------------| +| IntervalArithmetic | 0.22.12 | 1 | 1.0.4 | +| IntervalBoxes | 0.2 | 0.3 | 0.3.0 | +| IntervalContractors | 0.5 | 0.6 | 0.6.0 | +| ReversePropagation | 0.3 | 0.4 | 0.4.0 | +| StaticArrays | 1 | 1 | 1.9.18 (unchanged) | +| Symbolics | 5, 6 | 7 | 7.17.0 | + +### Code changes needed for compatibility + +**Removed `@register_symbolic x ∈ y::Interval`** from `src/IntervalConstraintProgramming.jl`: +- IntervalArithmetic v1.0 follows IEEE 1788 and deliberately does NOT define `Base.==` or `Base.isequal`/`Base.hash` for `Interval`. Users should use `isequal_interval` etc. +- SymbolicUtils uses `isequal`/`hash` for hash-consing symbolic expressions. Embedding an `Interval` in a symbolic expression (via `@register_symbolic x ∈ y::Interval`) triggers `isequal` → `==` → `InconclusiveBooleanOperation` error. +- **Decision: Do NOT define `Base.isequal`/`Base.hash` for `Interval`** — that would be type piracy contradicting IntervalArithmetic's design. Instead, remove the `∈` registration and avoid putting `Interval` values inside symbolic expressions. + +**Replaced `@register_symbolic x ∈ y::Interval`** with decomposition in `src/IntervalConstraintProgramming.jl`: +- New: `Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y)))` +- This decomposes `x ∈ a..b` into `(x >= a) & (x <= b)` at the symbolic level, avoiding Interval values in the symbolic tree entirely. +- Users can still write `x^2 + y^2 ∈ interval(0, 1)` — it just gets decomposed into two comparison constraints combined with `&`. + +**Changed `Separator` constructor** in `src/contractor.jl`: +- Old: `Separator(ex, vars, constraint::Interval) = Separator(vars, ex ∈ constraint, constraint, ...)` +- New: `Separator(ex, vars, constraint::Interval) = Separator(vars, ex, constraint, ...)` +- The `ex` field no longer wraps in `∈ constraint` (the constraint is already stored separately in the `constraint` field). + +**Updated `show` for `AbstractSeparator`** to display constraint info when available (via `hasproperty` check), since `ex` no longer contains it. + +**Fixed pre-existing bug in `separator()` in `src/utils.jl`**: +- The `&` and `|` handlers used `∩`/`∪` (`Base.intersect`/`Base.union`) instead of `⊓`/`⊔` (from IntervalArithmetic.Symbols, defined for separators in `set_operations.jl`). +- This bug was never triggered before because tests didn't exercise the `separator()` path with `&`/`|`. It surfaced now because `∈` decomposes into `(expr >= lo) & (expr <= hi)`. + +### Known limitations + +- 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`. diff --git a/Project.toml b/Project.toml index 978104c..e9536c2 100644 --- a/Project.toml +++ b/Project.toml @@ -11,12 +11,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] -IntervalArithmetic = "0.22.12" -IntervalBoxes = "0.2" -IntervalContractors = "0.5" -ReversePropagation = "0.3" +IntervalArithmetic = "1" +IntervalBoxes = "0.3" +IntervalContractors = "0.6" +ReversePropagation = "0.4" StaticArrays = "1" -Symbolics = "5, 6" +Symbolics = "7" julia = "1" [extras] diff --git a/future.md b/future.md new file mode 100644 index 0000000..76ba180 --- /dev/null +++ b/future.md @@ -0,0 +1,18 @@ +# Future work + +## `@constraint` macro for chained comparisons + +Julia lowers `0 <= x^2+y^2 <= 1` to `(0 <= x^2+y^2) && (x^2+y^2 <= 1)`, and `&&` requires a `Bool`, so this fails with symbolic expressions. + +A `@constraint` macro could intercept the AST before lowering and convert chained comparisons into `&` (which works symbolically) or directly into the `∈` form: + +```julia +@constraint 0 <= x^2 + y^2 <= 1 # → x^2+y^2 ∈ interval(0, 1) +@constraint x^2 + y^2 <= 1 # → simple case, works as-is +``` + +The package already exports `@constraint` but it is not yet defined. + +## ReversePropagation method overwrite warnings + +ReversePropagation emits many "Method definition overwritten" warnings with Symbolics v7. Harmless but noisy — should be addressed upstream in ReversePropagation. diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl index cca9f62..12b02be 100644 --- a/src/IntervalConstraintProgramming.jl +++ b/src/IntervalConstraintProgramming.jl @@ -24,12 +24,16 @@ import Base: import IntervalArithmetic: mid, interval, emptyinterval, isinf, isinterior, hull, mince - @register_symbolic ¬(x) -@register_symbolic x ∈ y::Interval @register_symbolic x ∨ y @register_symbolic x ∧ y +# We cannot register `x ∈ y::Interval` as a symbolic operation because +# SymbolicUtils hash-consing requires isequal/hash, which Interval deliberately +# does not define (IEEE 1788). Instead, decompose into comparisons that the +# package already handles. +Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y))) + export # BasicContractor, # @contractor, diff --git a/src/contractor.jl b/src/contractor.jl index 4aefbf2..3181bbd 100644 --- a/src/contractor.jl +++ b/src/contractor.jl @@ -25,7 +25,13 @@ end # Base.show(io::IO, S::Separator) = print(io, "Separator($(S.ex) ∈ $(S.constraint), vars = $(join(S.vars, ", ")))") -Base.show(io::IO, S::AbstractSeparator) = print(io, "Separator($(S.ex), vars=$(join(S.vars, ", ")))") +function Base.show(io::IO, S::AbstractSeparator) + if hasproperty(S, :constraint) + print(io, "Separator($(S.ex) ∈ $(S.constraint), vars=$(join(S.vars, ", ")))") + else + print(io, "Separator($(S.ex), vars=$(join(S.vars, ", ")))") + end +end function Separator(orig_expr, vars) ex, constraint = analyse(orig_expr) @@ -33,7 +39,7 @@ function Separator(orig_expr, vars) return Separator(ex, vars, constraint) end -Separator(ex, vars, constraint::Interval) = Separator(vars, ex ∈ constraint, constraint, make_function(ex, vars), Contractor(ex, vars)) +Separator(ex, vars, constraint::Interval) = Separator(vars, ex, constraint, make_function(ex, vars), Contractor(ex, vars)) function separate_infinite_box(S::Separator, X::IntervalBox) # for an box that extends to infinity we cannot evaluate at a corner diff --git a/src/utils.jl b/src/utils.jl index abff343..64b76cd 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -22,19 +22,19 @@ function analyse(ex) if op ∈ (≤, <) constraint = interval(-Inf, 0) Num(lhs - rhs), constraint - + elseif op ∈ (≥, >) constraint = interval(0, +Inf) Num(lhs - rhs), constraint - + elseif op == (==) constraint = interval(0, 0) Num(lhs - rhs), constraint - + else return ex, interval(0, 0) # implicit 0 end - + end @@ -59,23 +59,23 @@ function separator(ex, vars) lhs, rhs = arguments(ex2) if op == & - return separator(lhs, vars) ∩ separator(rhs, vars) + return separator(lhs, vars) ⊓ separator(rhs, vars) elseif op == | - return separator(lhs, vars) ∪ separator(rhs, vars) - + return separator(lhs, vars) ⊔ separator(rhs, vars) + elseif op ∈ (≤, <) constraint = interval(-Inf, 0) Separator(Num(lhs - rhs), vars, constraint) - + elseif op ∈ (≥, >) constraint = interval(0, +Inf) Separator(Num(lhs - rhs), vars, constraint) - + elseif op == (==) constraint = interval(0, 0) Separator(Num(lhs - rhs), vars, constraint) - + else return Separator(ex, vars, interval(0, 0)) # implicit "== 0" end From e54814deb1bed21ff5be30de67a15d17175fb058 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Thu, 2 Apr 2026 03:48:31 -0500 Subject: [PATCH 2/8] Bump version to 0.15.0 Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e9536c2..4ddf9aa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "IntervalConstraintProgramming" uuid = "138f1668-1576-5ad7-91b9-7425abbf3153" -version = "0.14.0" +version = "0.15.0" [deps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" From 6bed4b81cb185195c805ce22633eeb119cc2d323 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Thu, 2 Apr 2026 13:52:22 -0500 Subject: [PATCH 3/8] Update CI actions, require Julia 1.11, remove obsolete REQUIRE Cherry-pick infrastructure changes from PR #220: - Update GitHub Actions versions (checkout v6, setup-julia v2, cache v3, codecov v6) - Test on Julia 1.11 instead of 1.10 - Set julia compat to 1.11 - Remove obsolete REQUIRE file (Pkg.jl era) Co-Authored-By: Claude Opus 4.6 (1M context) --- .github/workflows/CI.yml | 10 +++++----- Project.toml | 2 +- REQUIRE | 6 ------ 3 files changed, 6 insertions(+), 12 deletions(-) delete mode 100644 REQUIRE diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 3b794d8..21c1b8a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -11,7 +11,7 @@ jobs: fail-fast: false matrix: version: - - '1.10' + - '1.11' - 'nightly' os: - ubuntu-latest @@ -24,12 +24,12 @@ jobs: - os: macOS-latest arch: x86 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v3 env: cache-name: cache-artifacts with: @@ -42,6 +42,6 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v6 with: file: lcov.info diff --git a/Project.toml b/Project.toml index 4ddf9aa..3da5f52 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,7 @@ IntervalContractors = "0.6" ReversePropagation = "0.4" StaticArrays = "1" Symbolics = "7" -julia = "1" +julia = "1.11" [extras] Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index f4b7bcc..0000000 --- a/REQUIRE +++ /dev/null @@ -1,6 +0,0 @@ -julia 1.0 -ModelingToolkit -IntervalArithmetic 0.15 -IntervalRootFinding 0.4 -IntervalContractors 0.3 -MacroTools 0.4 From 1783bb0ac7a95c4ef9ba4e3d0b13adbcfd9733de Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Thu, 2 Apr 2026 13:52:54 -0500 Subject: [PATCH 4/8] Test on Julia 1.10, 1.12, and nightly; require Julia >= 1.10 Co-Authored-By: Claude Opus 4.6 (1M context) --- .github/workflows/CI.yml | 3 ++- Project.toml | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 21c1b8a..ffd2714 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -11,7 +11,8 @@ jobs: fail-fast: false matrix: version: - - '1.11' + - '1.10' + - '1.12' - 'nightly' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index 3da5f52..e3fc06b 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,7 @@ IntervalContractors = "0.6" ReversePropagation = "0.4" StaticArrays = "1" Symbolics = "7" -julia = "1.11" +julia = "1.10" [extras] Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" From 299a3fc03a7b3d3cd0e83efefaba9898406a280b Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Thu, 2 Apr 2026 13:53:35 -0500 Subject: [PATCH 5/8] Widen compat ranges to also support older dependency versions Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index e3fc06b..d64eda0 100644 --- a/Project.toml +++ b/Project.toml @@ -11,12 +11,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] -IntervalArithmetic = "1" -IntervalBoxes = "0.3" -IntervalContractors = "0.6" -ReversePropagation = "0.4" +IntervalArithmetic = "0.22.12 - 0.23, 1" +IntervalBoxes = "0.2 - 0.3" +IntervalContractors = "0.5 - 0.6" +ReversePropagation = "0.3 - 0.4" StaticArrays = "1" -Symbolics = "7" +Symbolics = "5 - 7" julia = "1.10" [extras] From d5e943ca926a1f7f2222b510ee64c6a00c62ebff Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Mon, 13 Apr 2026 19:21:51 -0500 Subject: [PATCH 6/8] Set up documentation site with DocumenterVitepress Switch from the old Documenter HTML backend to DocumenterVitepress for a modern VitePress-based documentation site. Add new pages explaining contractors/separators and the internal architecture, update index.md to the current API, add GitHub Actions workflow for doc deployment, and remove stale mkdocs.yml and Manifest.toml. Co-Authored-By: Claude Opus 4.6 (1M context) --- .github/workflows/Documenter.yml | 47 +++++ README.md | 4 + docs/Manifest.toml | 145 --------------- docs/Project.toml | 6 + docs/make.jl | 32 ++-- docs/mkdocs.yml | 24 --- docs/package.json | 14 ++ docs/src/architecture.md | 157 ++++++++++++++++ docs/src/contractors_and_separators.md | 247 +++++++++++++++++++++++++ docs/src/index.md | 218 ++++------------------ 10 files changed, 528 insertions(+), 366 deletions(-) create mode 100644 .github/workflows/Documenter.yml delete mode 100644 docs/Manifest.toml delete mode 100644 docs/mkdocs.yml create mode 100644 docs/package.json create mode 100644 docs/src/architecture.md create mode 100644 docs/src/contractors_and_separators.md diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..67b78e0 --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,47 @@ +name: Documenter + +on: + push: + branches: + - master + tags: ['*'] + pull_request: + + workflow_dispatch: + +permissions: + contents: write + pages: write + id-token: write + statuses: write + +concurrency: + group: pages + cancel-in-progress: false + +jobs: + build: + runs-on: ubuntu-latest + steps: + - name: Checkout + uses: actions/checkout@v6 + - name: Setup Julia + uses: julia-actions/setup-julia@v2 + - name: Load Julia packages from cache + id: julia-cache + uses: julia-actions/cache@v3 + - name: Build and deploy docs + uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + GKSwstype: "100" + JULIA_DEBUG: "Documenter" + - name: Save Julia depot cache on cancel or failure + id: julia-cache-save + if: cancelled() || failure() + uses: actions/cache/save@v5 + with: + path: | + ${{ steps.julia-cache.outputs.cache-paths }} + key: ${{ steps.julia-cache.outputs.cache-key }} diff --git a/README.md b/README.md index 342dd71..bd3558e 100644 --- a/README.md +++ b/README.md @@ -52,6 +52,10 @@ plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary") +## Architecture + +See [docs/src/architecture.md](docs/src/architecture.md) for a detailed description of the internal pipeline from symbolic expressions to contractors and separators. + ## Author - [David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders), diff --git a/docs/Manifest.toml b/docs/Manifest.toml deleted file mode 100644 index 9f9b3b2..0000000 --- a/docs/Manifest.toml +++ /dev/null @@ -1,145 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -[[Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[CRlibm]] -deps = ["Libdl"] -git-tree-sha1 = "9d1c22cff9c04207f336b8e64840d0bd40d86e0e" -uuid = "96374032-68de-5a5b-8d9e-752f78720389" -version = "0.8.0" - -[[Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" - -[[DocStringExtensions]] -deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.1" - -[[Documenter]] -deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "646ebc3db49889ffeb4c36f89e5d82c6a26295ff" -uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.24.7" - -[[ErrorfreeArithmetic]] -git-tree-sha1 = "a5198ab6c8a724dd3965b31ddd11ccde65300f5b" -uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a" -version = "0.5.0" - -[[FastRounding]] -deps = ["ErrorfreeArithmetic", "Test"] -git-tree-sha1 = "224175e213fd4fe112db3eea05d66b308dc2bf6b" -uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f" -version = "0.2.0" - -[[InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[IntervalArithmetic]] -deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "SetRounding", "StaticArrays"] -git-tree-sha1 = "b2db6ee367b4eb3ee8b009ede8ca809e4fd23d35" -uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" -version = "0.16.7" - -[[JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" -uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.0" - -[[LibGit2]] -deps = ["Printf"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[LinearAlgebra]] -deps = ["Libdl"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[Parsers]] -deps = ["Dates", "Test"] -git-tree-sha1 = "75d07cb840c300084634b4991761886d0d762724" -uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.1" - -[[Pkg]] -deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" - -[[Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[Random]] -deps = ["Serialization"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[RecipesBase]] -git-tree-sha1 = "b4ed4a7f988ea2340017916f7c9e5d7560b52cae" -uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -version = "0.8.0" - -[[SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" - -[[Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[SetRounding]] -deps = ["Test"] -git-tree-sha1 = "faca28c7937138d560ede5bfef2076d0cdf3584b" -uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0" -version = "0.2.0" - -[[Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[SparseArrays]] -deps = ["LinearAlgebra", "Random"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - -[[StaticArrays]] -deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "5a3bcb6233adabde68ebc97be66e95dcb787424c" -uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "0.12.1" - -[[Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" - -[[Test]] -deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/docs/Project.toml b/docs/Project.toml index f64f5b3..d91698c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,9 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" + +[compat] +Documenter = "1" +DocumenterVitepress = "0.3" diff --git a/docs/make.jl b/docs/make.jl index 18f9c73..4e5f98e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,22 +1,28 @@ using Documenter +using DocumenterVitepress using IntervalConstraintProgramming, IntervalArithmetic -makedocs( +makedocs(; modules = [IntervalConstraintProgramming], - doctest = true, - format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"), authors = "David P. Sanders", sitename = "IntervalConstraintProgramming.jl", - - pages = Any[ + format = DocumenterVitepress.MarkdownVitepress(; + repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl", + devbranch = "master", + devurl = "dev", + ), + pages = [ "Home" => "index.md", - "API" => "api.md" - ] - ) + "Contractors and Separators" => "contractors_and_separators.md", + "Architecture" => "architecture.md", + "API" => "api.md", + ], +) -deploydocs( - repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl.git", - target = "build", - deps = nothing, - make = nothing, +DocumenterVitepress.deploydocs(; + repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl", + target = joinpath(@__DIR__, "build"), + branch = "gh-pages", + devbranch = "master", + push_preview = true, ) diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml deleted file mode 100644 index 860e60b..0000000 --- a/docs/mkdocs.yml +++ /dev/null @@ -1,24 +0,0 @@ -site_name: PACKAGE_NAME.jl -repo_url: https://github.com/USER_NAME/PACKAGE_NAME.jl -site_description: Description... -site_author: USER_NAME - -theme: readthedocs - -extra_css: - - assets/Documenter.css - -extra_javascript: - - https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML - - assets/mathjaxhelper.js - -markdown_extensions: - - extra - - tables - - fenced_code - - mdx_math - -docs_dir: 'build' - -pages: - - Home: index.md diff --git a/docs/package.json b/docs/package.json new file mode 100644 index 0000000..42eda9b --- /dev/null +++ b/docs/package.json @@ -0,0 +1,14 @@ +{ + "devDependencies": { + "@types/node": "^22.0.0" + }, + "scripts": { + "docs:dev": "vitepress dev build/.documenter", + "docs:build": "vitepress build build/.documenter", + "docs:preview": "vitepress preview build/.documenter" + }, + "dependencies": { + "vitepress": "^1.6.4", + "@mdit/plugin-mathjax": "^0.26.1" + } +} diff --git a/docs/src/architecture.md b/docs/src/architecture.md new file mode 100644 index 0000000..8147561 --- /dev/null +++ b/docs/src/architecture.md @@ -0,0 +1,157 @@ +# Architecture: From Symbolic Expressions to Contractors + +This document describes the internal pipeline that converts symbolic mathematical expressions into interval contractors and separators for constraint propagation. + +## Overview + +IntervalConstraintProgramming works in three stages: + +1. **Parse**: A symbolic constraint (e.g. `x^2 + y^2 <= 1`) is normalized to the form `f(x,y) in [a, b]`. +2. **Compile**: The expression `f` is decomposed into SSA (single static assignment) form via common subexpression elimination (CSE), then a forward--backward contractor is generated using [ReversePropagation.jl](https://github.com/JuliaIntervals/ReversePropagation.jl). +3. **Evaluate**: The separator applies the contractor to interval boxes, classifying regions as *inner* (definitely feasible), *outer* (definitely infeasible), or *boundary* (unknown). + +``` +Symbolic expression (Symbolics.jl) + | + separator(ex, vars) + | + v + Normalize: f(x) in constraint + | + +---> make_function(f, vars) --> evaluation function + +---> Contractor(f, vars) + | + v + forward_backward_contractor(f, vars) [ReversePropagation] + | + +---> cse_equations(f) --> SSA assignments + +---> constraint intersection + +---> rev.() on each assignment --> backward propagation + +---> eval() --> compiled Julia closure + | + v + Separator(vars, f, constraint, eval_fn, contractor) +``` + +## Stage 1: Parsing and normalization + +The entry point is `separator(ex)` (or equivalently `constraint(ex)`), defined in `src/utils.jl`. + +Given a symbolic expression, it: + +- **Detects the operator** (`<=`, `>=`, `==`, `&`, `|`, `!`) at the top level. +- **Handles logical connectives** recursively: + - `A & B` becomes `separator(A)` meet `separator(B)` + - `A | B` becomes `separator(A)` join `separator(B)` + - `!A` swaps inner and outer of `separator(A)` +- **Normalizes inequalities** by moving everything to one side: + - `lhs <= rhs` becomes `lhs - rhs` with constraint `[-Inf, 0]` + - `lhs >= rhs` becomes `lhs - rhs` with constraint `[0, +Inf]` + - `lhs == rhs` becomes `lhs - rhs` with constraint `[0, 0]` + +The `in` operator is handled symbolically: +```julia +Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y))) +``` +so `x^2 + y^2 in interval(0, 1)` decomposes into two inequality constraints combined with `&`. + +## Stage 2: Compilation via ReversePropagation + +### Common Subexpression Elimination (CSE) + +The function `cse_equations(ex)` in ReversePropagation (`src/cse.jl`) recursively decomposes the expression into binary operations, assigning a fresh variable to each subexpression: + +``` +Input: 3x^2 + 4x^2*y + +Output (SSA assignments): + _a := x * x + _b := 3 * _a + _c := _a * y + _d := 4 * _c + _e := _b + _d +``` + +Common subexpressions (like `x * x`) are computed only once. + +### Forward--backward contractor (HC4Revise) + +The function `forward_backward_code(ex, vars)` in `src/reverse_icp.jl` generates the full contractor: + +1. **Forward pass**: The CSE assignments evaluate `f` with interval arithmetic. +2. **Constraint intersection**: The result is intersected with the constraint interval: + ``` + value := _e + _e := _e meet constraint + ``` +3. **Backward pass**: Each CSE assignment is reversed using `rev()`, which calls the corresponding reverse function from [IntervalContractors.jl](https://github.com/JuliaIntervals/IntervalContractors.jl) (e.g. `plus_rev`, `mul_rev`, `sin_rev`). These propagate the constraint back through each operation, contracting the input intervals. + +The resulting code is compiled into a Julia closure via `eval()`. + +### Reverse operations + +For each elementary operation, IntervalContractors provides a *reverse* (or *backward*) function: + +| Forward | Reverse | Effect | +|-----------------|-----------------|-------------------------------------------| +| `z := x + y` | `plus_rev` | Contract `x` and `y` given constraint on `z` | +| `z := x * y` | `mul_rev` | Contract `x` and `y` given constraint on `z` | +| `z := sin(x)` | `sin_rev` | Contract `x` given constraint on `z` | +| `z := x ^ n` | `power_rev` | Contract `x` given constraint on `z` | + +These are registered as symbolic functions so they can appear in generated code. + +## Stage 3: Separator evaluation + +A `Separator` stores both the compiled evaluation function and the contractor. When called with an `IntervalBox`: + +```julia +function (SS::Separator)(X::IntervalBox) + # 1. Contract: find the boundary region + boundary = SS.contractor(X) + + # 2. Evaluate at corners to classify inner/outer + lb_image = SS.f(inf.(X)) + ub_image = SS.f(sup.(X)) + + # 3. If a corner satisfies the constraint, extend inner; otherwise extend outer + inner = boundary + outer = boundary + # ... extend with corner evaluations ... + + return (boundary, inner, outer) +end +``` + +For boxes extending to infinity, a different strategy is used: the contractor is applied once with the constraint and once with its complement, and inner/outer are determined from those two contractions. + +### Separator composition + +Separators compose via set-theoretic operations defined in `src/set_operations.jl`: + +| Operation | Symbol | Inner | Outer | +|-----------|--------|----------------------|----------------------| +| Intersect | `S1 meet S2` | `inner1 meet inner2` | `outer1 join outer2` | +| Union | `S1 join S2` | `inner1 join inner2` | `outer1 meet outer2` | +| Complement| `!S` | swap inner/outer | swap inner/outer | + +## Paving + +The `pave(X, S, epsilon)` function in `src/pave.jl` recursively subdivides the domain box `X`: + +1. Apply the separator to get `(boundary, inner, outer)`. +2. Regions proven inner are accumulated. +3. Regions proven outer are discarded. +4. Boundary regions are bisected and re-processed until their diameter falls below `epsilon`. + +The result is two collections: `inner` boxes (guaranteed feasible) and `boundary` boxes (undetermined at this tolerance). + +## Gradient computation + +ReversePropagation also provides `gradient_code(ex, vars)` which generates SSA code for the gradient of `ex` using reverse-mode automatic differentiation (via ChainRules.jl): + +1. **Forward pass**: Same CSE decomposition as above. +2. **Initialization**: Set the adjoint of the output to 1. +3. **Adjoint pass**: For each forward assignment in reverse order, compute the adjoint contribution using `rrule` from ChainRules. + +The result is a sequence of assignments whose final variables hold the partial derivatives. This SSA code can itself be fed back into the forward--backward machinery to create contractors for derivative constraints. diff --git a/docs/src/contractors_and_separators.md b/docs/src/contractors_and_separators.md new file mode 100644 index 0000000..1153749 --- /dev/null +++ b/docs/src/contractors_and_separators.md @@ -0,0 +1,247 @@ +# Contractors and Separators + +This page explains the key concepts behind IntervalConstraintProgramming.jl: **contractors** and **separators**. These are the tools that allow us to rigorously determine which regions of space satisfy a given set of constraints. + +## Interval boxes + +An **interval box** (or simply "box") is a Cartesian product of intervals, representing a rectangular region in $n$-dimensional space. For example, in 2D: + +$$X = [1, 3] \times [0, 2]$$ + +is the rectangle with $x \in [1, 3]$ and $y \in [0, 2]$. + +In Julia, we write: + +```julia +using IntervalArithmetic, IntervalBoxes +X = IntervalBox(interval(1, 3), interval(0, 2)) +``` + +## Contractors + +A **contractor** $C$ is an operator that takes an interval box $X$ and returns a smaller (or equal) box $C(X) \subseteq X$. The key property is: + +> **Contractance**: $C(X) \subseteq X$ for all boxes $X$. + +Given an expression like $x^2 + y^2$ and a constraint interval like $[0, 1]$, a contractor shrinks the input box to eliminate regions that are *guaranteed* not to satisfy the constraint $x^2 + y^2 \in [0, 1]$. + +### How contraction works: forward--backward propagation + +IntervalConstraintProgramming uses a **forward--backward contractor**, also known as constraint propagation: + +1. **Forward pass**: Evaluate the expression using interval arithmetic. Each intermediate result gets an interval that encloses all possible values. + +2. **Constraint intersection**: Intersect the result with the constraint interval. If the expression evaluates to $[{-0.5}, 1.5]$ and the constraint is $[0, 1]$, the intersected output is $[0, 1]$. + +3. **Backward pass**: Propagate the tightened output *backwards* through each operation, using inverse interval operations. This tightens the intervals on the input variables. + +For example, if we know $z = x + y$, $z \in [0, 1]$, $x \in [{-1}, 3]$ and $y \in [{-1}, 3]$, then the backward step deduces $x \in [{-1}, 2]$ and $y \in [{-1}, 2]$, since $x$ cannot exceed $1 - ({-1}) = 2$ given the constraint on $z$. + +### Creating a contractor + +```julia +using IntervalArithmetic, IntervalConstraintProgramming, Symbolics + +vars = @variables x, y +C = Contractor(x^2 + y^2, vars) +``` + +Calling the contractor on a box returns the contracted box: + +```julia +X = IntervalBox(-5..5, 2) + +C(X) # contract with default constraint [0, 0] +C(X, interval(0, 1)) # contract with constraint [0, 1] +``` + +Contractors are useful on their own when you want to find the set where an expression takes a particular value (or range of values). However, they do not distinguish the *inside* of a set from the *outside* --- that is what separators do. + + +## Separators + +A **separator** $S$ extends the idea of a contractor to classify regions as **inner** (satisfying the constraint), **outer** (violating the constraint), or **boundary** (undetermined). + +Given a constraint like $x^2 + y^2 \le 1$, a separator applied to a box $X$ returns three boxes: + +- **inner**: a box guaranteed to lie *inside* the feasible set (where the constraint is satisfied) +- **outer**: a box guaranteed to contain the *outside* (where the constraint is violated) +- **boundary**: a box containing the boundary of the constraint set, where the status is unknown + +``` +boundary, inner, outer = S(X) +``` + +The fundamental guarantee is: + +> Every point in $X$ that satisfies the constraint lies in `inner`. +> Every point in $X$ that violates the constraint lies in `outer`. +> The `boundary` box contains all points whose status is uncertain. + +### Creating a separator + +The main way to create a separator is with the `constraint` function (or equivalently `separator`): + +```julia +using IntervalArithmetic, IntervalArithmetic.Symbols +using IntervalConstraintProgramming, Symbolics + +vars = @variables x, y + +S = constraint(x^2 + y^2 <= 1, vars) +``` + +The constraint is automatically normalized: `x^2 + y^2 <= 1` becomes `x^2 + y^2 - 1` with constraint interval $[-\infty, 0]$. + +You can use any of the standard comparison operators: + +| Constraint syntax | Meaning | +|---------------------------------|-------------------------------------------| +| `constraint(expr <= val, vars)` | $\text{expr} \le \text{val}$ | +| `constraint(expr >= val, vars)` | $\text{expr} \ge \text{val}$ | +| `constraint(expr ~ val, vars)` | $\text{expr} = \text{val}$ (equality) | + +You can also use the `in` operator ($\in$) with an interval: + +```julia +S = constraint(x^2 + y^2 in interval(1, 3), vars) +``` + +This is equivalent to $1 \le x^2 + y^2 \le 3$. + +### Applying a separator + +```julia +X = IntervalBox(-5..5, 2) + +boundary, inner, outer = S(X) +``` + +For a large box like $[-5, 5]^2$, the contractor can only tighten slightly. The real power comes from combining separators with the **paving algorithm**. + +## Paving: finding the feasible set + +The `pave` function recursively bisects the domain and applies a separator at each step. Boxes that are proved to be inside the constraint set are collected as the **inner** approximation; boxes that cannot be classified at the given tolerance are collected as the **boundary**: + +```julia +X = IntervalBox(-5..5, 2) +S = constraint(x^2 + y^2 <= 1, vars) + +inner, boundary = pave(X, S, 0.05) +``` + +The third argument is the **tolerance**: bisection stops when a boundary box has diameter smaller than this value. + +The result gives rigorous guarantees: +- Every point in an `inner` box satisfies the constraint. +- Every point *not* covered by `inner` or `boundary` boxes violates the constraint. +- Decreasing the tolerance gives a sharper approximation at the cost of more computation. + +### Plotting the result + +```julia +using Plots + +plot(collect.(inner), aspectratio=1, lw=0, label="inner") +plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary") +``` + +## Combining separators: set operations + +Separators can be combined using set-theoretic operations to describe more complex feasible sets. + +### Intersection ($\sqcap$) + +The intersection of two separators gives the set of points satisfying *both* constraints: + +```julia +S1 = constraint(x^2 + y^2 >= 1, vars) +S2 = constraint(x^2 + y^2 <= 3, vars) + +S = S1 ⊓ S2 # type \sqcap in Julia +``` + +This defines the annular region $1 \le x^2 + y^2 \le 3$. + +The inner contractor of `S1 ⊓ S2` is the intersection of the inner contractors, and the outer contractor is the union of the outer contractors: + +$$\text{inner}(S_1 \sqcap S_2) = \text{inner}(S_1) \sqcap \text{inner}(S_2)$$ +$$\text{outer}(S_1 \sqcap S_2) = \text{outer}(S_1) \sqcup \text{outer}(S_2)$$ + +### Union ($\sqcup$) + +The union gives the set of points satisfying *either* constraint: + +```julia +S = S1 ⊔ S2 # type \sqcup in Julia +``` + +$$\text{inner}(S_1 \sqcup S_2) = \text{inner}(S_1) \sqcup \text{inner}(S_2)$$ +$$\text{outer}(S_1 \sqcup S_2) = \text{outer}(S_1) \sqcap \text{outer}(S_2)$$ + + +### Complement (`!`) + +The complement swaps inner and outer: + +```julia +S_complement = !S +``` + +### Set difference and symmetric difference + +```julia +S_diff = setdiff(S1, S2) # S1 and not S2 +S_sym = symdiff(S1, S2) # in one but not both +``` + +These are defined in terms of the basic operations: +- `setdiff(S1, S2)` = `S1 ⊓ !S2` +- `symdiff(S1, S2)` = `setdiff(S1, S2) ⊔ setdiff(S2, S1)` + +## Complete example + +Here is a complete example that finds the intersection of an ellipse and a rotated constraint: + +```julia +using IntervalArithmetic, IntervalArithmetic.Symbols +using IntervalConstraintProgramming +using IntervalBoxes +using Symbolics + +vars = @variables x, y + +# An elliptical region +C1 = constraint(x^2 + 2y^2 >= 1, vars) + +# A second constraint +C2 = constraint(x^2 + y^2 + x * y <= 3, vars) + +# Intersection: points satisfying both constraints +C = C1 ⊓ C2 + +# Pave the domain +X = IntervalBox(-5..5, 2) +inner, boundary = pave(X, C, 0.05) + +# Visualize +using Plots +plot(collect.(inner), aspectratio=1, lw=0, label="inner") +plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary") +``` + +## How it all fits together + +The pipeline from a symbolic constraint to a paving is: + +1. **Symbolic expression** $\to$ normalized form $f(\mathbf{x}) \in [a, b]$ +2. **Forward--backward contractor** compiled from the expression via [ReversePropagation.jl](https://github.com/JuliaIntervals/ReversePropagation.jl) +3. **Separator** wraps the contractor and classifies regions as inner, outer, or boundary +4. **Paving algorithm** recursively bisects and applies the separator to build a rigorous approximation of the feasible set + +For more details on the internal compilation pipeline, see the [Architecture](architecture.md) page. + +## References + +- Jaulin, L., Kieffer, M., Didrit, O., and Walter, E. (2001). *Applied Interval Analysis*. Springer. +- Jaulin, L. and Desrochers, B. (2014). Introduction to the Algebra of Separators with Application to Path Planning. *Engineering Applications of Artificial Intelligence*, **33**, 141--147. diff --git a/docs/src/index.md b/docs/src/index.md index 3dd5ede..d2d4686 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,214 +1,64 @@ -# `IntervalConstraintProgramming.jl` +# IntervalConstraintProgramming.jl This Julia package allows you to specify a set of **constraints** on real-valued variables, given by (in)equalities, and rigorously calculate inner and outer approximations to the *feasible set*, i.e. the set that satisfies the constraints. -This uses interval arithmetic provided by the author's -[`IntervalArithmetic.jl`](https://github.com/dpsanders/IntervalArithmetic.jl) package, +This uses interval arithmetic provided by +[`IntervalArithmetic.jl`](https://github.com/JuliaIntervals/IntervalArithmetic.jl), in particular multi-dimensional `IntervalBox`es, i.e. Cartesian products of one-dimensional intervals. To do this, *interval constraint programming* is used, in particular the so-called "forward--backward contractor". This is implemented in terms of *separators*; see -[Jaulin & Desrochers]. +[Jaulin & Desrochers (2014)](https://doi.org/10.1016/j.engappai.2014.04.010). -```@meta -DocTestSetup = quote - using IntervalConstraintProgramming, IntervalArithmetic -end -``` - -## Usage -Let's define a constraint, using the `@constraint` macro: -```jldoctest -julia> using IntervalConstraintProgramming, IntervalArithmetic - -julia> S = @constraint x^2 + y^2 <= 1 -Separator: - - variables: x, y - - expression: x ^ 2 + y ^ 2 ∈ [-Inf, 1] -``` -It works out automatically that `x` and `y` are variables. -The macro creates a `Separator` object, in this case a `ConstraintSeparator`. - -We now create an initial interval box in the ``x``--``y`` plane: -```julia -julia> x = y = -100..100 # notation for creating an interval with `IntervalArithmetic.jl` - -julia> X = IntervalBox(x, y) -``` - -The `@constraint` macro defines an object `S`, of type `Separator`. -This is a function which, -when applied to the box ``X = x \times y`` -in the x--y plane, applies two *contractors*, an inner one and an outer one. - -The inner contractor tries to shrink down, or *contract*, the box, to the smallest subbox -of ``X`` that contains the part of ``X`` that satisfies the constraint; the -outer contractor tries to contract ``X`` to the smallest subbox that contains the -region where the constraint is not satisfied. - -When `S` is applied to the box `X`, it returns the result of the inner and outer contractors: -```julia -julia> inner, outer = S(X); - -julia> inner -([-1, 1],[-1, 1]) - -julia> outer -([-100, 100],[-100, 100]) -``` - -### Without using macros - -We can also make an object `S`, of type `Separator` or `C`, of type `Contractor` without using Macros, for that you need to define variables using `ModelingToolkit.jl`. -Example - -```julia -julia> using IntervalConstraintProgramming, IntervalArithmetic, ModelingToolkit - -julia> @variables x y -(x(), y()) - -julia> S = Separator(x+y<1) -Separator: - - variables: x, y - - expression: x() + y() == [-Inf, 1] - -julia> C = Contractor(x+y) - Contractor in 2 dimensions: - - forward pass contracts to 1 dimensions - - variables: Symbol[:x, :y] - - expression: x() + y() -``` - -While making `Separator`or `Contractor`'s object we can also specify variables, like this - -```julia -julia> vars = @variables x y z -(x(), y(), z()) - -julia> S = Separator(vars, x+y<1) -Separator: - - variables: x, y, z - - expression: x() + y() == [-Inf, 1] - -julia> C = Contractor(vars, y+z) -Contractor in 3 dimensions: - - forward pass contracts to 1 dimensions - - variables: Symbol[:x, :y, :z] - - expression: y() + z() -``` -We can make objects (of type `Separator` or `Contractor`)by just using function name (Note: you have to specify variables explicitly as discussed above when you make objects by using function name). We can also use polynomial function to make objects. +## Quick start ```julia -julia> vars=@variables x y -(x(), y()) +using IntervalArithmetic, IntervalArithmetic.Symbols +using IntervalConstraintProgramming +using IntervalBoxes +using Symbolics -julia> f(a,b)= a+b -f (generic function with 1 method) +vars = @variables x, y -julia> C = Contractor(vars,f) -Contractor in 2 dimensions: - - forward pass contracts to 1 dimensions - - variables: Symbol[:x, :y] - - expression: x() + y() +# Create separators from constraints +C1 = constraint(x^2 + 2y^2 >= 1, vars) +C2 = constraint(x^2 + y^2 + x * y <= 3, vars) -julia> f(a,b) = a+b <1 -f (generic function with 1 method) +# Combine with intersection +C = C1 ⊓ C2 -julia> S=Separator(vars, f) -Separator: - - variables: x, y - - expression: x() + y() == [-Inf, 1] - -julia> using DynamicPolynomials #using polynomial functions - -julia> pvars = @polyvar x y -(x, y) - -julia> f(a,b) = a + b < 1 -p (generic function with 1 method) - -julia> S=Separator(pvars, f) -Separator: - - variables: x, y - - expression: x() + y() == [-Inf, 1] -``` -#### BasicContractor -Objects of type `Contractor` have four fields (variables, forward, backward and expression), among them data of two fields (forward, backward) are useful (i.e forward and backward functions) for further usage of that object, thats why it is preferred to use an object of type `BasicContractor` in place of `Contractor` which only contain these two fields for less usage of memory by unloading all the extra stuff.(Note: Like object of `Contractor` type,`BasicContractor`'s object will also have all the properties which are discussed above). - -```julia -julia> @variables x y -(x(), y()) - -julia> C = BasicContractor(x^2 + y^2) - Basic version of Contractor +# Pave the domain +X = IntervalBox(-5..5, 2) +inner, boundary = pave(X, C, 0.05) ``` - - -## Set inversion: finding the feasible set - -To make progress, we must recursively bisect and apply the contractors, keeping -track of the region proved to be inside the feasible set, and the region that is -on the boundary ("both inside and outside"). This is done by the `pave` function, -that takes a separator, a domain to search inside, and an optional tolerance: +## Plotting the result ```julia -julia> using Plots +using Plots -julia> x = y = -100..100 - -julia> S = @constraint 1 <= x^2 + y^2 <= 3 - -julia> paving = pave(S, X, 0.125); +plot(collect.(inner), aspectratio=1, lw=0, label="inner") +plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary") ``` -`pave` returns an object of type `Paving`. This contains: the separator itself; -an `inner` approximation, of type `SubPaving`, which is an alias for a `Vector` of `IntervalBox`es; -a `SubPaving` representing the boxes on the boundary that could not be assigned either to the inside or outside of the set; -and the tolerance. - -We may draw the result using a plot recipe from `IntervalArithmetic`. Either a -single `IntervalBox`, or a `Vector` of `IntervalBox`es (which a `SubPaving` is) -maybe be drawn using `plot` from `Plots.jl`: -```julia -julia> plot(paving.inner, c="green") -julia> plot!(paving.boundary, c="gray") -``` - -The output should look something like this: - -![Ring](ring.png) - - -The green boxes have been **rigorously** proved to be contained within the feasible set, -and the white boxes to be outside the set. The grey boxes show those that lie on the boundary, whose status is unknown. - -### 3D - -The package works in any number of dimensions, although it suffers from the usual exponential slowdown in higher dimensions ("combinatorial explosion"); in 3D, it is still relatively fast. - -There are sample 3D calculations in the `examples` directory, in particular in the [solid torus notebook](examples/Solid torus.ipynb), which uses [`GLVisualize.gl`](https://github.com/JuliaGL/GLVisualize.jl) to provide an interactive visualization that may be rotated and zoomed. The output for the solid torus looks like this: - -![Coloured solid torus](solid_torus.png) - +The inner (blue) region is **rigorously** guaranteed to lie *inside* the feasible set. +The white region is guaranteed to lie *outside* the set. +The boundary (red) region lies on the boundary, where the status is unknown at the given tolerance. ## Set operations -Separators may be combined using the operators `!` (complement), `⊓` and `⊔` to make -more complicated sets; see the [notebook](https://github.com/JuliaIntervals/IntervalConstraintProgrammingNotebooks/blob/master/Basic%20examples%20of%20separators.ipynb) for several examples. Further examples can be found in the repository [IntervalConstraintProgrammingNotebooks](https://github.com/JuliaIntervals/IntervalConstraintProgrammingNotebooks). -## Author +Separators may be combined using the operators `!` (complement), `⊓` (intersection) and `⊔` (union) to build +more complicated sets. See the [Contractors and Separators](@ref) page for details. -- [David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders) - - [Julia lab, MIT](http://julia.mit.edu/) - - Departamento de Física, Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM) +## Learn more +- [Contractors and Separators](contractors_and_separators.md) -- detailed explanation of the key concepts +- [Architecture](architecture.md) -- how the internal pipeline works +- [API](api.md) -- function and type reference -## References: -- *Applied Interval Analysis*, Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter (2001) -- Introduction to the Algebra of Separators with Application to Path Planning, Luc Jaulin and Benoît Desrochers, *Engineering Applications of Artificial Intelligence* **33**, 141–147 (2014) +## References -## Acknowledements -Financial support is acknowledged from DGAPA-UNAM PAPIME grants PE-105911 and PE-107114, and DGAPA-UNAM PAPIIT grant IN-117214, and from a CONACYT-Mexico sabbatical fellowship. The author thanks Alan Edelman and the Julia group for hospitality during his sabbatical visit. He also thanks Luc Jaulin and Jordan Ninin for the [IAMOOC](http://iamooc.ensta-bretagne.fr/) online course, which introduced him to this subject. +- *Applied Interval Analysis*, Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter (2001) +- Introduction to the Algebra of Separators with Application to Path Planning, Luc Jaulin and Benoît Desrochers, *Engineering Applications of Artificial Intelligence* **33**, 141--147 (2014) From d9b615f117db0d0c339691edb4771e6a70253a7b Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Tue, 14 Apr 2026 13:25:16 -0500 Subject: [PATCH 7/8] Add DAG-based constraint propagation as alternative to code generation Implements an explicit DAG (directed acyclic graph) for forward-backward interval constraint propagation (HC4Revise), inspired by Schichl & Neumaier. This provides an inspectable, iterable alternative to the existing ReversePropagation code-generation approach. New types: DAGContractor, DAGSeparator, ConstraintDAG New files: src/dag/{nodes,build,propagate,contractor}.jl Tests: test/test_dag.jl (27 tests) Benchmarks: benchmark/bench_dag_vs_codegen.jl Both approaches produce identical paving results. The DAG approach is currently ~10x slower due to DAG reconstruction per call and dynamic dispatch, with clear optimization paths documented in CLAUDE.md. Co-Authored-By: Claude Opus 4.6 (1M context) --- CLAUDE.md | 56 +++++++ benchmark/bench_dag_vs_codegen.jl | 76 +++++++++ src/IntervalConstraintProgramming.jl | 9 + src/dag/build.jl | 99 +++++++++++ src/dag/contractor.jl | 119 +++++++++++++ src/dag/nodes.jl | 67 ++++++++ src/dag/propagate.jl | 239 +++++++++++++++++++++++++++ test/test_dag.jl | 149 +++++++++++++++++ 8 files changed, 814 insertions(+) create mode 100644 benchmark/bench_dag_vs_codegen.jl create mode 100644 src/dag/build.jl create mode 100644 src/dag/contractor.jl create mode 100644 src/dag/nodes.jl create mode 100644 src/dag/propagate.jl create mode 100644 test/test_dag.jl diff --git a/CLAUDE.md b/CLAUDE.md index 395e1fc..a219bc1 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -43,3 +43,59 @@ 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`. + +## DAG-based constraint propagation (2026-04-14) + +### Overview + +Added an alternative constraint propagation engine based on explicit DAGs (directed acyclic graphs), inspired by the Schichl & Neumaier approach. Lives in `src/dag/` alongside the existing code-generation approach. + +**New files:** +- `src/dag/nodes.jl` — `AbstractNode`, `VariableNode`, `ConstantNode`, `OperationNode`, `ConstraintDAG` +- `src/dag/build.jl` — `build_dag()`: walks Symbolics expression tree, builds DAG with CSE via `objectid`-keyed cache. Decomposes n-ary `+`/`*` into left-associated binary chains. +- `src/dag/propagate.jl` — `forward!`, `backward!`, `propagate!` (iterative HC4Revise) +- `src/dag/contractor.jl` — `DAGContractor`, `DAGSeparator` types +- `test/test_dag.jl` — 27 tests (construction, propagation, pave, comparison with original) +- `benchmark/bench_dag_vs_codegen.jl` — head-to-head benchmarks + +**New exports:** `DAGContractor`, `DAGSeparator`, `ConstraintDAG` + +### Architecture decisions + +**DAG rebuilt per contractor call.** The `DAGContractor` rebuilds the DAG each time it's called with a (possibly different) constraint interval. This is correct but slower than caching. The `Separator` protocol requires contracting with `f(x) = 0` (boundary) by default, plus `f(x) ∈ [a,b]` and complement for infinite boxes — so the same expression needs different constraints. A future optimization would cache DAGs per constraint value. + +**Mutable node ranges.** Each node stores a mutable `range::Interval` that is updated in-place during propagation. This avoids allocating new interval containers per pass. The trade-off is that the DAG is not thread-safe. + +**Reverse contractor dispatch.** Uses explicit `_contract_binary!`/`_contract_unary!` methods dispatching on `typeof(op)` (e.g., `::typeof(+)`, `::typeof(sin)`). Unary ops are generated via `@eval` loop over the IntervalContractors reverse functions. Falls back with a warning for unknown ops. + +**N-ary decomposition.** Symbolics.jl represents `x + y + 1` as `+(x, y, 1)` with 3 arguments. The DAG builder decomposes these into left-associated binary chains: `(x + y) + 1`. This is necessary because `plus_rev` and `mul_rev` are binary. + +### Performance (benchmarks) + +The DAG approach produces **identical results** (same box counts) but is ~10x slower due to: +1. DAG reconstruction per call (no caching) +2. Dynamic dispatch on `node.op` for each forward/backward step +3. Allocation of child-range vectors in `_eval_forward` + +| Problem | ϵ | Codegen | DAG | Ratio | +|---------|---|---------|-----|-------| +| Unit disk 2D | 0.1 | 287 μs | 3.2 ms | 11x | +| Unit disk 2D | 0.01 | 2.3 ms | 27 ms | 12x | +| Annulus 2D | 0.1 | 1.6 ms | 14 ms | 9x | +| 3D torus | 1.0 | 12 ms | 120 ms | 10x | + +### Potential optimizations (not yet implemented) + +1. **Cache DAGs** by constraint value in the separator (avoid rebuilding each call) +2. **Typed operation enum** instead of `Function` field to eliminate dynamic dispatch +3. **Pre-allocated workspace** for child ranges (avoid per-node allocation) +4. **Compiled propagation schedule** — a flat array of `(op_code, input_indices, output_index)` instructions that a tight loop can execute without dynamic dispatch, similar to what the codegen approach produces but without `eval()` +5. **Multi-constraint propagation** — multiple constraints sharing variable nodes in a single DAG, so narrowing from one constraint immediately benefits others (the main theoretical advantage of the Schichl-Neumaier approach over independent contractors) + +### Findings during implementation + +**Separator boundary semantics.** The existing `Separator` calls its contractor with `constraint = interval(0.0)` (contracting onto the boundary `f(x) = 0`), not with the full constraint interval. This is how it distinguishes inner from boundary from outer. Initial DAG implementation incorrectly contracted with the full constraint, producing no inner boxes. Fixed by making `DAGContractor` accept the constraint as a call-time argument (defaulting to `interval(0.0)`). + +**Symbolics literal representation.** Numbers in Symbolics expression trees (e.g., the `2` in `x^2`, the `3` in `3*x`) are wrapped as `BasicSymbolic{SymReal}` objects, not Julia `Number`s. They satisfy `!issym && !iscall` and can be unwrapped via `Symbolics.value()` to get the underlying `Int64`/`Float64`. Variable names are extracted via `Symbolics.nameof()`. + +**`+` and `*` parse issue.** Writing `op === + || op === *` causes a Julia parse error because `||` tries to parse `*` as a boolean expression. Must parenthesize: `op === (+) || op === (*)`. diff --git a/benchmark/bench_dag_vs_codegen.jl b/benchmark/bench_dag_vs_codegen.jl new file mode 100644 index 0000000..19a5695 --- /dev/null +++ b/benchmark/bench_dag_vs_codegen.jl @@ -0,0 +1,76 @@ +using IntervalConstraintProgramming +using IntervalArithmetic, IntervalArithmetic.Symbols +using IntervalBoxes +using Symbolics +using BenchmarkTools + +@variables x, y, z + +println("=" ^ 60) +println("Benchmark: DAG vs code-generation constraint propagation") +println("=" ^ 60) + +# --- Problem 1: Unit disk (x² + y² ≤ 1) --- + +println("\n--- Problem 1: Unit disk (x² + y² ≤ 1) ---\n") + +S_codegen = Separator(x^2 + y^2 <= 1, [x, y]) +S_dag = DAGSeparator(x^2 + y^2 <= 1, [x, y]) +X = IntervalBox(-3..3, 2) + +println("Single separator call on [-3,3]²:") +print(" codegen: ") +@btime $S_codegen($X) +print(" DAG: ") +@btime $S_dag($X) + +for ϵ in [1.0, 0.1, 0.01] + println("\npave with ϵ = $ϵ:") + print(" codegen: ") + r1 = @btime pave($X, $S_codegen, $ϵ) + print(" DAG: ") + r2 = @btime pave($X, $S_dag, $ϵ) + inner1, bnd1 = r1 + inner2, bnd2 = r2 + println(" codegen: $(length(inner1)) inner, $(length(bnd1)) boundary") + println(" DAG: $(length(inner2)) inner, $(length(bnd2)) boundary") +end + +# --- Problem 2: Annulus (1 ≤ x² + y² ≤ 4) --- + +println("\n--- Problem 2: Annulus (1 ≤ x² + y² ≤ 4) ---\n") + +annulus_expr = x^2 + y^2 ∈ interval(1, 4) +S_codegen2 = separator(annulus_expr, [x, y]) +S_dag2 = DAGSeparator(annulus_expr, [x, y]) + +X2 = IntervalBox(-5..5, 2) +for ϵ in [1.0, 0.1] + println("pave with ϵ = $ϵ:") + print(" codegen: ") + r1 = @btime pave($X2, $S_codegen2, $ϵ) + print(" DAG: ") + r2 = @btime pave($X2, $S_dag2, $ϵ) + inner1, bnd1 = r1 + inner2, bnd2 = r2 + println(" codegen: $(length(inner1)) inner, $(length(bnd1)) boundary") + println(" DAG: $(length(inner2)) inner, $(length(bnd2)) boundary") +end + +# --- Problem 3: 3D constraint --- + +println("\n--- Problem 3: 3D (3 - sqrt(x²+y²))² + z² ≤ 1) ---\n") + +S_codegen3 = Separator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z]) +S_dag3 = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z]) + +X3 = IntervalBox(-10..10, 3) +println("pave with ϵ = 1.0:") +print(" codegen: ") +r1 = @btime pave($X3, $S_codegen3, 1.0) +print(" DAG: ") +r2 = @btime pave($X3, $S_dag3, 1.0) +inner1, bnd1 = r1 +inner2, bnd2 = r2 +println(" codegen: $(length(inner1)) inner, $(length(bnd1)) boundary") +println(" DAG: $(length(inner2)) inner, $(length(bnd2)) boundary") diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl index 12b02be..0b46af8 100644 --- a/src/IntervalConstraintProgramming.jl +++ b/src/IntervalConstraintProgramming.jl @@ -54,6 +54,9 @@ export export ¬ +# DAG-based propagation types +export DAGContractor, DAGSeparator, ConstraintDAG + const reverse_operations = IntervalContractors.reverse_operations @@ -62,5 +65,11 @@ include("contractor.jl") include("set_operations.jl") include("pave.jl") +# DAG-based alternative to code-generation approach +include("dag/nodes.jl") +include("dag/build.jl") +include("dag/propagate.jl") +include("dag/contractor.jl") + end # module diff --git a/src/dag/build.jl b/src/dag/build.jl new file mode 100644 index 0000000..48379be --- /dev/null +++ b/src/dag/build.jl @@ -0,0 +1,99 @@ +const SU = Symbolics.SymbolicUtils + +""" + build_dag(expr::Num, vars::Vector{Num}, constraint::Interval) → ConstraintDAG + +Build a DAG from a Symbolics expression. Uses a cache dictionary for +common subexpression elimination (CSE): identical symbolic subexpressions +map to the same node. + +N-ary `+` and `*` are decomposed into left-associated binary operations +so that each node has exactly 2 children for binary ops (or 1 for unary). +""" +function build_dag(expr::Num, variables::Vector{Num}, constraint::Interval) + unwrapped_vars = Symbolics.value.(variables) + sym_expr = Symbolics.value(expr) + + # Map variable symbols to VariableNodes + var_nodes = VariableNode[] + var_map = Dict{Any, VariableNode}() + for (i, v) in enumerate(unwrapped_vars) + vn = VariableNode(Symbolics.nameof(v), i, entireinterval()) + push!(var_nodes, vn) + var_map[v] = vn + end + + # Cache for CSE: symbolic subexpr → AbstractNode + cache = Dict{UInt, AbstractNode}() + all_nodes = AbstractNode[] + + root = _build_node!(sym_expr, var_map, cache, all_nodes) + + return ConstraintDAG(var_nodes, all_nodes, root, constraint) +end + +""" +Recursively build DAG nodes from a symbolic expression. +Returns the node corresponding to `expr`. +""" +function _build_node!(expr, var_map, cache, all_nodes) + # Check cache first (CSE) + h = objectid(expr) + haskey(cache, h) && return cache[h] + + node = if haskey(var_map, expr) + # It's a known variable + vn = var_map[expr] + if vn ∉ all_nodes + push!(all_nodes, vn) + end + vn + + elseif Symbolics.issym(expr) + # Unknown symbol — shouldn't happen if vars are correctly specified + error("Unknown variable in expression: $(Symbolics.nameof(expr))") + + elseif SU.iscall(expr) + _build_operation!(expr, var_map, cache, all_nodes) + + else + # Literal number + val = Symbolics.value(expr) + cn = ConstantNode(val) + push!(all_nodes, cn) + cn + end + + cache[h] = node + return node +end + +""" +Build an OperationNode from a symbolic function call. +Handles n-ary `+` and `*` by decomposing to left-associated binary chains. +""" +function _build_operation!(expr, var_map, cache, all_nodes) + op = Symbolics.operation(expr) + args = Symbolics.arguments(expr) + + if (op === (+) || op === (*)) && length(args) > 2 + # Decompose n-ary to binary: (((a + b) + c) + d) + node = _build_node!(args[1], var_map, cache, all_nodes) + for i in 2:length(args) + right = _build_node!(args[i], var_map, cache, all_nodes) + parent = OperationNode(op, AbstractNode[node, right]) + push!(all_nodes, parent) + node = parent + end + return node + end + + # Standard case: unary or binary operation + child_nodes = AbstractNode[ + _build_node!(a, var_map, cache, all_nodes) for a in args + ] + + node = OperationNode(op, child_nodes) + push!(all_nodes, node) + return node +end diff --git a/src/dag/contractor.jl b/src/dag/contractor.jl new file mode 100644 index 0000000..6ab9a40 --- /dev/null +++ b/src/dag/contractor.jl @@ -0,0 +1,119 @@ +""" + DAGContractor <: AbstractContractor + +A contractor based on explicit DAG propagation (HC4Revise). +Alternative to the code-generation approach in `Contractor`. + +Stores the expression and variables so it can build a DAG with any constraint +at call time. By default contracts onto `f(x) = 0` (the boundary), matching +the behavior of the code-generation `Contractor`. +""" +struct DAGContractor{V} <: AbstractContractor + vars::V + ex::Any +end + +DAGContractor(ex::Num, vars::Vector{Num}) = DAGContractor(vars, ex) + +function (C::DAGContractor)(X::IntervalBox, constraint::Interval=interval(0.0)) + dag = build_dag(C.ex, C.vars, constraint) + return propagate!(dag, X) +end + +""" + DAGSeparator <: AbstractSeparator + +A separator that uses DAG-based propagation. +Drop-in replacement for `Separator` — same `(boundary, inner, outer)` output. +""" +struct DAGSeparator{V, E, C, F} <: AbstractSeparator + vars::V + ex::E + constraint::C + f::F + contractor::DAGContractor +end + +function DAGSeparator(expr::Num, vars::Vector{Num}, constraint::Interval) + contractor = DAGContractor(expr, vars) + f = make_function(expr, vars) + return DAGSeparator(vars, expr, constraint, f, contractor) +end + +function (SS::DAGSeparator)(X) + if any(x -> isinf(diam(x)), X) + return _separate_infinite_dag(SS, X) + end + + # Contract onto the boundary (f(x) = 0), same as original Separator + boundary = SS.contractor(X) + + lb = IntervalBox(inf.(X)) + ub = IntervalBox(sup.(X)) + + inner = boundary + outer = boundary + + lb_image = SS.f(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) + if !isempty_interval(ub_image) && issubset_interval(ub_image, SS.constraint) + inner = inner ⊔ ub + else + outer = outer ⊔ ub + end + + return boundary, inner, outer +end + +function _separate_infinite_dag(SS::DAGSeparator, X::IntervalBox) + C = SS.contractor + a, b = inf(SS.constraint), sup(SS.constraint) + + inner = C(X, interval(a, b)) + + local outer + if a == -Inf + outer = C(X, interval(b, Inf)) + elseif b == Inf + outer = C(X, interval(-Inf, a)) + else + outer1 = C(X, interval(-Inf, a)) + outer2 = C(X, interval(b, Inf)) + outer = outer1 ⊔ outer2 + end + + boundary = inner ⊓ outer + return (boundary, inner, outer) +end + +""" + DAGSeparator(orig_expr, vars) + +Build a DAG-based separator from a symbolic constraint expression, +handling boolean combinations (`&`, `|`, `!`) by composing separators. +""" +function DAGSeparator(orig_expr, vars) + ex_val = Symbolics.value(orig_expr) + + if SU.iscall(ex_val) + op = Symbolics.operation(ex_val) + args = Symbolics.arguments(ex_val) + + if op === (&) + return DAGSeparator(Num(args[1]), vars) ⊓ DAGSeparator(Num(args[2]), vars) + elseif op === (|) + return DAGSeparator(Num(args[1]), vars) ⊔ DAGSeparator(Num(args[2]), vars) + elseif op === (!) + return ¬(DAGSeparator(Num(args[1]), vars)) + end + end + + ex, constraint = analyse(orig_expr) + return DAGSeparator(ex, vars, constraint) +end diff --git a/src/dag/nodes.jl b/src/dag/nodes.jl new file mode 100644 index 0000000..8f577c4 --- /dev/null +++ b/src/dag/nodes.jl @@ -0,0 +1,67 @@ +abstract type AbstractNode end + +mutable struct VariableNode <: AbstractNode + name::Symbol + index::Int # position in the input IntervalBox + range::Interval # current interval bound +end + +mutable struct ConstantNode <: AbstractNode + value::Float64 + range::Interval # [value, value], set once +end + +ConstantNode(v::Real) = ConstantNode(Float64(v), interval(Float64(v))) + +mutable struct OperationNode <: AbstractNode + op::Function # +, *, sin, ^, etc. + children::Vector{AbstractNode} + range::Interval # current interval bound (set during forward pass) +end + +OperationNode(op, children) = OperationNode(op, children, entireinterval()) + +# Convenience accessors + +is_variable(n::VariableNode) = true +is_variable(::AbstractNode) = false + +is_constant(n::ConstantNode) = true +is_constant(::AbstractNode) = false + +is_operation(n::OperationNode) = true +is_operation(::AbstractNode) = false + +function Base.show(io::IO, n::VariableNode) + print(io, "Var($(n.name), $(n.range))") +end + +function Base.show(io::IO, n::ConstantNode) + print(io, "Const($(n.value))") +end + +function Base.show(io::IO, n::OperationNode) + print(io, "Op($(n.op), $(length(n.children)) children, $(n.range))") +end + +""" + ConstraintDAG + +Directed acyclic graph representing an expression with a constraint on its output. + +- `variables`: input variable nodes (ordered to match IntervalBox) +- `nodes`: all nodes in topological order (leaves first, root last) +- `root`: the output node +- `constraint`: interval the root must lie in +""" +struct ConstraintDAG + variables::Vector{VariableNode} + nodes::Vector{AbstractNode} # topological order + root::AbstractNode + constraint::Interval +end + +function Base.show(io::IO, dag::ConstraintDAG) + vars = join([string(v.name) for v in dag.variables], ", ") + print(io, "ConstraintDAG($(length(dag.nodes)) nodes, vars=($vars), constraint=$(dag.constraint))") +end diff --git a/src/dag/propagate.jl b/src/dag/propagate.jl new file mode 100644 index 0000000..50511bd --- /dev/null +++ b/src/dag/propagate.jl @@ -0,0 +1,239 @@ +""" + forward!(dag::ConstraintDAG, X::IntervalBox) + +Forward propagation: evaluate intervals from leaves to root in topological order. +Sets `range` on every node. +""" +function forward!(dag::ConstraintDAG, X::IntervalBox) + for v in dag.variables + v.range = X[v.index] + end + + for node in dag.nodes + if node isa OperationNode + node.range = _eval_forward(node.op, node.children) + end + # VariableNode and ConstantNode ranges are already set + end + + return dag.root.range +end + +function _eval_forward(op::Function, children::Vector{AbstractNode}) + if length(children) == 1 + return op(children[1].range) + elseif length(children) == 2 + return op(children[1].range, children[2].range) + else + # Shouldn't happen after n-ary decomposition, but handle gracefully + return op([c.range for c in children]...) + end +end + +# Special case: ^ with integer exponent. Symbolics stores ^(x, 2) where +# the second arg is a ConstantNode. IntervalArithmetic dispatches x^n. +function _eval_forward(::typeof(^), children::Vector{AbstractNode}) + base = children[1].range + exp_node = children[2] + if exp_node isa ConstantNode + n = exp_node.value + if n == round(Int, n) + return base ^ round(Int, n) + end + end + return base ^ exp_node.range +end + + +""" + backward!(dag::ConstraintDAG) + +Backward propagation: propagate constraint from root to leaves in reverse +topological order using reverse interval contractors. + +The root's range is first intersected with the constraint. Then each +OperationNode propagates narrowed ranges to its children. +""" +function backward!(dag::ConstraintDAG) + dag.root.range = IntervalArithmetic.intersect_interval(dag.root.range, dag.constraint) + + if isempty_interval(dag.root.range) + # Constraint is infeasible on this box — mark everything empty + for v in dag.variables + v.range = emptyinterval() + end + return + end + + # Traverse in reverse topological order (root → leaves) + for i in length(dag.nodes):-1:1 + node = dag.nodes[i] + node isa OperationNode || continue + _contract_backward!(node) + end +end + +""" +Apply the reverse contractor for a single operation node, narrowing +both `node.range` and each child's range in place. +""" +function _contract_backward!(node::OperationNode) + op = node.op + children = node.children + z = node.range + + if isempty_interval(z) + for c in children + c.range = emptyinterval() + end + return + end + + if length(children) == 1 + _contract_unary!(op, z, children[1]) + elseif length(children) == 2 + _contract_binary!(op, z, children[1], children[2]) + end +end + +# --- Binary reverse contractors --- + +function _contract_binary!(::typeof(+), z, left, right) + z_new, l_new, r_new = plus_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +function _contract_binary!(::typeof(-), z, left, right) + z_new, l_new, r_new = minus_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +function _contract_binary!(::typeof(*), z, left, right) + z_new, l_new, r_new = mul_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +function _contract_binary!(::typeof(/), z, left, right) + z_new, l_new, r_new = div_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +function _contract_binary!(::typeof(^), z, base_node, exp_node) + if exp_node isa ConstantNode + n = exp_node.value + if n == round(Int, n) + # power_rev(z, x, n::Integer) → (z', x', n) + z_new, x_new, _ = power_rev(z, base_node.range, round(Int, n)) + base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new) + return + end + end + # General case: x^y + z_new, x_new, y_new = power_rev(z, base_node.range, exp_node.range) + base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new) + exp_node.range = IntervalArithmetic.intersect_interval(exp_node.range, y_new) +end + +function _contract_binary!(::typeof(min), z, left, right) + z_new, l_new, r_new = min_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +function _contract_binary!(::typeof(max), z, left, right) + z_new, l_new, r_new = max_rev(z, left.range, right.range) + left.range = IntervalArithmetic.intersect_interval(left.range, l_new) + right.range = IntervalArithmetic.intersect_interval(right.range, r_new) +end + +# --- Unary reverse contractors --- + +for (jl_op, rev_fn) in [ + (:sin, :sin_rev), + (:cos, :cos_rev), + (:tan, :tan_rev), + (:asin, :asin_rev), + (:acos, :acos_rev), + (:atan, :atan_rev), + (:sinh, :sinh_rev), + (:cosh, :cosh_rev), + (:tanh, :tanh_rev), + (:asinh, :asinh_rev), + (:acosh, :acosh_rev), + (:atanh, :atanh_rev), + (:exp, :exp_rev), + (:exp2, :exp2_rev), + (:exp10, :exp10_rev), + (:expm1, :expm1_rev), + (:log, :log_rev), + (:log2, :log2_rev), + (:log10, :log10_rev), + (:log1p, :log1p_rev), + (:sqrt, :sqrt_rev), + (:abs, :abs_rev), + (:sign, :sign_rev), +] + @eval function _contract_unary!(::typeof($jl_op), z, child) + z_new, x_new = $rev_fn(z, child.range) + child.range = IntervalArithmetic.intersect_interval(child.range, x_new) + end +end + +# Unary minus: -(x) +function _contract_unary!(::typeof(-), z, child) + z_new, x_new = minus_rev(z, child.range) + child.range = IntervalArithmetic.intersect_interval(child.range, x_new) +end + +# Fallback for unknown operations — no contraction, just leave ranges as-is +function _contract_unary!(op, z, child) + @warn "No reverse contractor for unary $op — skipping" maxlog=1 +end + +function _contract_binary!(op, z, left, right) + @warn "No reverse contractor for binary $op — skipping" maxlog=1 +end + + +""" + propagate!(dag::ConstraintDAG, X::IntervalBox; max_iter=10, tol=1e-10) → IntervalBox + +Iterate forward-backward propagation on the DAG until convergence or `max_iter` +iterations. Returns the narrowed IntervalBox. + +This implements the HC4Revise algorithm: each iteration does one forward pass +(natural interval evaluation) followed by one backward pass (constraint narrowing). +For expressions where a variable appears multiple times, a single pass cannot +fully propagate — iteration is needed because each pass tightens the shared node, +which then feeds tighter bounds into the next forward evaluation. +""" +function propagate!(dag::ConstraintDAG, X::IntervalBox; max_iter::Int=10, tol::Float64=1e-10) + for _ in 1:max_iter + old_ranges = IntervalBox(Tuple(v.range for v in dag.variables)) + + forward!(dag, X) + backward!(dag) + + # Build new X from narrowed variable ranges + X = IntervalBox(Tuple(v.range for v in dag.variables)) + + if isempty(X) + return X + end + + # Check convergence: did the variable ranges change significantly? + max_change = maximum( + diam(old) - diam(new) + for (old, new) in zip(old_ranges, X) + ) + if max_change < tol + break + end + end + + return X +end diff --git a/test/test_dag.jl b/test/test_dag.jl new file mode 100644 index 0000000..cbf62ee --- /dev/null +++ b/test/test_dag.jl @@ -0,0 +1,149 @@ +using IntervalArithmetic, IntervalArithmetic.Symbols +using IntervalConstraintProgramming +using IntervalConstraintProgramming: build_dag, forward!, backward!, propagate!, ConstraintDAG +using Symbolics +using IntervalBoxes +using Test + +const IntervalType{T} = Union{Interval{T}, BareInterval{T}} +eq(a::IntervalType, b::IntervalType) = isequal_interval(bareinterval(a), bareinterval(b)) +eq(a::IntervalBox, b::IntervalBox) = all(eq.(a, b)) +eq(a::Vector, b::Vector) = length(a) == length(b) && all(eq.(a, b)) + +@testset "DAG construction" begin + vars = @variables x, y + + dag = build_dag(x + y, vars, interval(0.0)) + @test length(dag.variables) == 2 + @test dag.variables[1].name == :x + @test dag.variables[2].name == :y + + dag2 = build_dag(x^2 + y^2, vars, interval(-Inf, 0.0)) + # Should have: x, y, 2 (const), x^2, y^2, x^2+y^2 = 6 nodes min + @test length(dag2.nodes) >= 6 + + # N-ary: x + y + 1 should decompose to binary + dag3 = build_dag(x + y + 1, vars, interval(0.0)) + ops = filter(n -> n isa IntervalConstraintProgramming.OperationNode, dag3.nodes) + # Each + node should have exactly 2 children + for op in ops + @test length(op.children) <= 2 + end +end + +@testset "DAG forward propagation" begin + vars = @variables x, y + + dag = build_dag(x^2 + y^2, vars, interval(-Inf, 1.0)) + X = IntervalBox(interval(0.0, 1.0), interval(0.0, 1.0)) + result = forward!(dag, X) + # x^2 ∈ [0, 1], y^2 ∈ [0, 1], sum ∈ [0, 2] but with the -1 from analyse... + # Actually analyse(x^2+y^2 <= 1) gives expr = -1 + x^2 + y^2, constraint = (-∞, 0] + # Let's test the raw expression directly + @test inf(result) <= 2.0 + @test sup(result) >= 0.0 +end + +@testset "DAG propagation contracts correctly" begin + vars = @variables x, y + + # x^2 + y^2 - 1 ≤ 0 (the unit disk) + ex, constraint = IntervalConstraintProgramming.analyse(x^2 + y^2 <= 1) + dag = build_dag(ex, vars, constraint) + X = IntervalBox(interval(-10.0, 10.0), 2) + result = propagate!(dag, X) + # Should narrow to [-1, 1]² + @test all(inf.(result) .>= -1.0 - 1e-10) + @test all(sup.(result) .<= 1.0 + 1e-10) + + # Box already inside: should not widen + X2 = IntervalBox(interval(0.0, 0.3), 2) + dag2 = build_dag(ex, vars, constraint) + result2 = propagate!(dag2, X2) + @test all(inf.(result2) .>= 0.0 - 1e-10) + @test all(sup.(result2) .<= 0.3 + 1e-10) +end + +@testset "DAGContractor" begin + vars = @variables x, y + + C = DAGContractor(x^2 + y^2, vars) + X = IntervalBox(-Inf..Inf, 2) + result = C(X, -Inf..1) + @test eq(result, IntervalBox(-1..1, 2)) + + X2 = IntervalBox(interval(0.5, 1.5), 3) + vars3 = @variables x, y, z + C1 = DAGContractor(x + y, vars3) + @test eq(C1(X2, -Inf..1), IntervalBox(interval(0.5, 0.5), interval(0.5, 0.5), interval(0.5, 1.5))) + + C2 = DAGContractor(y + z, vars3) + @test eq(C2(X2, -Inf..1), IntervalBox(interval(0.5, 1.5), interval(0.5, 0.5), interval(0.5, 0.5))) +end + +@testset "DAGSeparator" begin + vars = @variables x, y + + II = -100..100 + X = IntervalBox(II, 2) + S = DAGSeparator(x^2 + y^2 <= 1, vars) + + boundary, inner, outer = S(X) + @test eq(inner, IntervalBox(-1..1, 2)) + @test eq(outer, IntervalBox(II, 2)) + + X_inf = IntervalBox(-Inf..Inf, -Inf..Inf) + boundary2, inner2, outer2 = S(X_inf) + @test eq(inner2, IntervalBox(-1..1, 2)) + @test eq(outer2, IntervalBox(-Inf..Inf, 2)) +end + +@testset "DAG pave" begin + vars = @variables x, y + + S = DAGSeparator(x^2 + y^2 <= 1, vars) + X = IntervalBox(-3..3, 2) + inner, boundary = pave(X, S, 0.5) + @test length(inner) > 0 + @test length(boundary) > 0 + + inner_area = sum(prod(diam.(b)) for b in inner) + total_area = inner_area + sum(prod(diam.(b)) for b in boundary) + @test inner_area > 2.0 # should be close to π + @test total_area < 5.0 # should bracket π from above +end + +@testset "DAG pave matches original" begin + vars = @variables x, y + + S_orig = Separator(x^2 + y^2 <= 1, vars) + S_dag = DAGSeparator(x^2 + y^2 <= 1, vars) + + X = IntervalBox(-3..3, 2) + inner_orig, boundary_orig = pave(X, S_orig, 0.5) + inner_dag, boundary_dag = pave(X, S_dag, 0.5) + + @test length(inner_orig) == length(inner_dag) + @test length(boundary_orig) == length(boundary_dag) +end + +@testset "DAG with combined separators" begin + vars = @variables x, y + + S1 = DAGSeparator(x > 0, vars) + S2 = DAGSeparator(y > 0, vars) + S = S1 ⊓ S2 + + X = IntervalBox(-3..3, 2) + boundary, inner, outer = S(X) + @test eq(inner, IntervalBox(0..3, 2)) +end + +@testset "DAG 3D torus" begin + vars = @variables x, y, z + + S = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, vars) + X = IntervalBox(-10..10, 3) + inner, boundary = pave(X, S, 1.0) + @test inner[1] isa IntervalBox{3, Float64, Interval{Float64}} +end From 662b57c71806f76c4e8f40ee99e9b2e87ca10c35 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Tue, 14 Apr 2026 16:11:06 -0500 Subject: [PATCH 8/8] Make DAG persistent with shared nodes across constraints Redesign the DAG engine so the graph is built once and reused: - SharedDAG: persistent DAG that accumulates constraints via add_expression!. Variable nodes and common subexpressions (CSE) are shared across all expressions added to the same DAG. - Multi-constraint propagation: forward-backward passes contract all constraints jointly, so narrowing from one constraint immediately benefits others through shared variable nodes. - DAGContractor/DAGSeparator accept an existing SharedDAG, allowing multiple separators to share the same graph. Benchmarks show ~30% speedup vs the previous rebuild-per-call approach, and shared DAGs use 33% fewer nodes than separate ones. Resolves Project.toml merge conflict from cherry-pick. Co-Authored-By: Claude Opus 4.6 (1M context) --- CLAUDE.md | 63 ++++------ benchmark/bench_dag_vs_codegen.jl | 39 +++++- src/IntervalConstraintProgramming.jl | 2 +- src/dag/build.jl | 94 ++++++--------- src/dag/contractor.jl | 51 +++++--- src/dag/nodes.jl | 55 ++++++--- src/dag/propagate.jl | 174 +++++++++++++-------------- test/test_dag.jl | 126 +++++++++++-------- 8 files changed, 335 insertions(+), 269 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index a219bc1..c7c257d 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -48,54 +48,37 @@ This file tracks decisions and context discovered while working on this package. ### Overview -Added an alternative constraint propagation engine based on explicit DAGs (directed acyclic graphs), inspired by the Schichl & Neumaier approach. Lives in `src/dag/` alongside the existing code-generation approach. +Added an alternative constraint propagation engine based on explicit, persistent, shared DAGs (directed acyclic graphs), inspired by the Schichl & Neumaier approach. Lives in `src/dag/`. -**New files:** -- `src/dag/nodes.jl` — `AbstractNode`, `VariableNode`, `ConstantNode`, `OperationNode`, `ConstraintDAG` -- `src/dag/build.jl` — `build_dag()`: walks Symbolics expression tree, builds DAG with CSE via `objectid`-keyed cache. Decomposes n-ary `+`/`*` into left-associated binary chains. -- `src/dag/propagate.jl` — `forward!`, `backward!`, `propagate!` (iterative HC4Revise) -- `src/dag/contractor.jl` — `DAGContractor`, `DAGSeparator` types -- `test/test_dag.jl` — 27 tests (construction, propagation, pave, comparison with original) -- `benchmark/bench_dag_vs_codegen.jl` — head-to-head benchmarks +**Key types:** +- `SharedDAG` — persistent DAG accumulating constraints. Multiple expressions share variable nodes and common subexpressions (CSE via `objectid`-keyed cache). Built incrementally with `add_expression!`. +- `DAGContractor` / `DAGSeparator` — drop-in replacements for `Contractor`/`Separator` backed by a SharedDAG. Can share a DAG (`DAGContractor(dag, expr, vars)`) or create their own. +- `ConstraintEntry` — a `(root_node, constraint_interval)` pair stored in the SharedDAG. -**New exports:** `DAGContractor`, `DAGSeparator`, `ConstraintDAG` +**Files:** `src/dag/{nodes,build,propagate,contractor}.jl`, `test/test_dag.jl`, `benchmark/bench_dag_vs_codegen.jl` -### Architecture decisions +### Architecture -**DAG rebuilt per contractor call.** The `DAGContractor` rebuilds the DAG each time it's called with a (possibly different) constraint interval. This is correct but slower than caching. The `Separator` protocol requires contracting with `f(x) = 0` (boundary) by default, plus `f(x) ∈ [a,b]` and complement for infinite boxes — so the same expression needs different constraints. A future optimization would cache DAGs per constraint value. +The DAG is **persistent**: built once, reused across all `pave` iterations. The same `SharedDAG` can back multiple contractors/separators. When `add_expression!` is called, existing nodes (variables and common subexpressions like `x^2`) are reused rather than duplicated. -**Mutable node ranges.** Each node stores a mutable `range::Interval` that is updated in-place during propagation. This avoids allocating new interval containers per pass. The trade-off is that the DAG is not thread-safe. +Propagation with multiple constraints: `propagate!(dag, X)` does forward evaluation (leaves→root), then backward contraction (root→leaves) for **all** constraints in the DAG. Narrowing from one constraint immediately benefits the others since they share variable nodes. -**Reverse contractor dispatch.** Uses explicit `_contract_binary!`/`_contract_unary!` methods dispatching on `typeof(op)` (e.g., `::typeof(+)`, `::typeof(sin)`). Unary ops are generated via `@eval` loop over the IntervalContractors reverse functions. Falls back with a warning for unknown ops. +For separator boundary/complement contraction, `propagate!(dag, X, constraint)` overrides the stored constraint on the first root, allowing the same DAG to be used for boundary (`f(x)=0`), inner (`f(x)∈[a,b]`), and complement contraction without rebuilding. -**N-ary decomposition.** Symbolics.jl represents `x + y + 1` as `+(x, y, 1)` with 3 arguments. The DAG builder decomposes these into left-associated binary chains: `(x + y) + 1`. This is necessary because `plus_rev` and `mul_rev` are binary. +### Performance (benchmarks, persistent shared DAG) -### Performance (benchmarks) +| Problem | ϵ | Codegen | DAG (persistent) | Ratio | +|---------|---|---------|-----------------|-------| +| Unit disk 2D | single call | 1.2 μs | 6.8 μs | 6x | +| Unit disk 2D | 0.1 | 294 μs | 2.5 ms | 8x | +| Unit disk 2D | 0.01 | 2.5 ms | 20 ms | 8x | +| Annulus 2D | 0.1 | 1.7 ms | 10 ms | 6x | +| 3D torus | 1.0 | 12 ms | 97 ms | 8x | -The DAG approach produces **identical results** (same box counts) but is ~10x slower due to: -1. DAG reconstruction per call (no caching) -2. Dynamic dispatch on `node.op` for each forward/backward step -3. Allocation of child-range vectors in `_eval_forward` +Multi-constraint shared DAG: 10 nodes shared vs 15 separate (33% reduction). Joint propagation of `x²+y²≤1` and `x²-y≤0` correctly narrows `[-10,10]²` to `[-1,1]×[0,1]`. -| Problem | ϵ | Codegen | DAG | Ratio | -|---------|---|---------|-----|-------| -| Unit disk 2D | 0.1 | 287 μs | 3.2 ms | 11x | -| Unit disk 2D | 0.01 | 2.3 ms | 27 ms | 12x | -| Annulus 2D | 0.1 | 1.6 ms | 14 ms | 9x | -| 3D torus | 1.0 | 12 ms | 120 ms | 10x | +### Remaining optimizations -### Potential optimizations (not yet implemented) - -1. **Cache DAGs** by constraint value in the separator (avoid rebuilding each call) -2. **Typed operation enum** instead of `Function` field to eliminate dynamic dispatch -3. **Pre-allocated workspace** for child ranges (avoid per-node allocation) -4. **Compiled propagation schedule** — a flat array of `(op_code, input_indices, output_index)` instructions that a tight loop can execute without dynamic dispatch, similar to what the codegen approach produces but without `eval()` -5. **Multi-constraint propagation** — multiple constraints sharing variable nodes in a single DAG, so narrowing from one constraint immediately benefits others (the main theoretical advantage of the Schichl-Neumaier approach over independent contractors) - -### Findings during implementation - -**Separator boundary semantics.** The existing `Separator` calls its contractor with `constraint = interval(0.0)` (contracting onto the boundary `f(x) = 0`), not with the full constraint interval. This is how it distinguishes inner from boundary from outer. Initial DAG implementation incorrectly contracted with the full constraint, producing no inner boxes. Fixed by making `DAGContractor` accept the constraint as a call-time argument (defaulting to `interval(0.0)`). - -**Symbolics literal representation.** Numbers in Symbolics expression trees (e.g., the `2` in `x^2`, the `3` in `3*x`) are wrapped as `BasicSymbolic{SymReal}` objects, not Julia `Number`s. They satisfy `!issym && !iscall` and can be unwrapped via `Symbolics.value()` to get the underlying `Int64`/`Float64`. Variable names are extracted via `Symbolics.nameof()`. - -**`+` and `*` parse issue.** Writing `op === + || op === *` causes a Julia parse error because `||` tries to parse `*` as a boolean expression. Must parenthesize: `op === (+) || op === (*)`. +1. **Typed operation enum** instead of `Function` field to eliminate dynamic dispatch +2. **Pre-allocated workspace** for child ranges +3. **Compiled propagation schedule** — flat instruction array for tight-loop execution diff --git a/benchmark/bench_dag_vs_codegen.jl b/benchmark/bench_dag_vs_codegen.jl index 19a5695..c805808 100644 --- a/benchmark/bench_dag_vs_codegen.jl +++ b/benchmark/bench_dag_vs_codegen.jl @@ -1,4 +1,5 @@ using IntervalConstraintProgramming +using IntervalConstraintProgramming: SharedDAG, add_expression!, propagate! using IntervalArithmetic, IntervalArithmetic.Symbols using IntervalBoxes using Symbolics @@ -7,7 +8,7 @@ using BenchmarkTools @variables x, y, z println("=" ^ 60) -println("Benchmark: DAG vs code-generation constraint propagation") +println("Benchmark: DAG (persistent) vs code-generation") println("=" ^ 60) # --- Problem 1: Unit disk (x² + y² ≤ 1) --- @@ -16,7 +17,7 @@ println("\n--- Problem 1: Unit disk (x² + y² ≤ 1) ---\n") S_codegen = Separator(x^2 + y^2 <= 1, [x, y]) S_dag = DAGSeparator(x^2 + y^2 <= 1, [x, y]) -X = IntervalBox(-3..3, 2) +X = IntervalBox(interval(-3, 3), 2) println("Single separator call on [-3,3]²:") print(" codegen: ") @@ -44,7 +45,7 @@ annulus_expr = x^2 + y^2 ∈ interval(1, 4) S_codegen2 = separator(annulus_expr, [x, y]) S_dag2 = DAGSeparator(annulus_expr, [x, y]) -X2 = IntervalBox(-5..5, 2) +X2 = IntervalBox(interval(-5, 5), 2) for ϵ in [1.0, 0.1] println("pave with ϵ = $ϵ:") print(" codegen: ") @@ -64,7 +65,7 @@ println("\n--- Problem 3: 3D (3 - sqrt(x²+y²))² + z² ≤ 1) ---\n") S_codegen3 = Separator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z]) S_dag3 = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z]) -X3 = IntervalBox(-10..10, 3) +X3 = IntervalBox(interval(-10, 10), 3) println("pave with ϵ = 1.0:") print(" codegen: ") r1 = @btime pave($X3, $S_codegen3, 1.0) @@ -74,3 +75,33 @@ inner1, bnd1 = r1 inner2, bnd2 = r2 println(" codegen: $(length(inner1)) inner, $(length(bnd1)) boundary") println(" DAG: $(length(inner2)) inner, $(length(bnd2)) boundary") + +# --- Problem 4: Multi-constraint on shared DAG --- + +println("\n--- Problem 4: Multi-constraint shared DAG ---") +println(" x²+y² ≤ 1 AND x²-y ≤ 0 (jointly propagated)\n") + +dag_shared = SharedDAG([x, y]) +add_expression!(dag_shared, x^2 + y^2 - 1, interval(-Inf, 0.0)) +add_expression!(dag_shared, x^2 - y, interval(-Inf, 0.0)) + +dag_sep1 = SharedDAG([x, y]) +add_expression!(dag_sep1, x^2 + y^2 - 1, interval(-Inf, 0.0)) +dag_sep2 = SharedDAG([x, y]) +add_expression!(dag_sep2, x^2 - y, interval(-Inf, 0.0)) + +X4 = IntervalBox(interval(-10, 10), 2) + +println("Single propagation call on [-10,10]²:") +print(" shared DAG (joint): ") +r_shared = @btime propagate!($dag_shared, $X4) +print(" separate DAGs (seq): ") +r_sep = @btime begin + r1 = propagate!($dag_sep1, $X4) + propagate!($dag_sep2, r1) +end +println(" shared result: ", r_shared) +println(" separate result: ", r_sep) + +println("\nShared DAG nodes: ", length(dag_shared.nodes), + " vs separate: ", length(dag_sep1.nodes) + length(dag_sep2.nodes)) diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl index 0b46af8..b4bad11 100644 --- a/src/IntervalConstraintProgramming.jl +++ b/src/IntervalConstraintProgramming.jl @@ -55,7 +55,7 @@ export export ¬ # DAG-based propagation types -export DAGContractor, DAGSeparator, ConstraintDAG +export DAGContractor, DAGSeparator, SharedDAG, add_expression! const reverse_operations = IntervalContractors.reverse_operations diff --git a/src/dag/build.jl b/src/dag/build.jl index 48379be..f756bd6 100644 --- a/src/dag/build.jl +++ b/src/dag/build.jl @@ -1,99 +1,83 @@ const SU = Symbolics.SymbolicUtils """ - build_dag(expr::Num, vars::Vector{Num}, constraint::Interval) → ConstraintDAG + add_expression!(dag::SharedDAG, expr::Num, constraint::Interval) → ConstraintEntry -Build a DAG from a Symbolics expression. Uses a cache dictionary for -common subexpression elimination (CSE): identical symbolic subexpressions -map to the same node. +Add an expression with its constraint to an existing SharedDAG. Reuses +variable nodes and any common subexpressions already in the DAG. -N-ary `+` and `*` are decomposed into left-associated binary operations -so that each node has exactly 2 children for binary ops (or 1 for unary). +Returns the ConstraintEntry that was added. """ -function build_dag(expr::Num, variables::Vector{Num}, constraint::Interval) - unwrapped_vars = Symbolics.value.(variables) +function add_expression!(dag::SharedDAG, expr::Num, constraint::Interval) sym_expr = Symbolics.value(expr) + root = _build_node!(sym_expr, dag) + entry = ConstraintEntry(root, constraint) + push!(dag.constraints, entry) + return entry +end - # Map variable symbols to VariableNodes - var_nodes = VariableNode[] - var_map = Dict{Any, VariableNode}() - for (i, v) in enumerate(unwrapped_vars) - vn = VariableNode(Symbolics.nameof(v), i, entireinterval()) - push!(var_nodes, vn) - var_map[v] = vn - end - - # Cache for CSE: symbolic subexpr → AbstractNode - cache = Dict{UInt, AbstractNode}() - all_nodes = AbstractNode[] - - root = _build_node!(sym_expr, var_map, cache, all_nodes) +""" + build_dag(expr::Num, vars::Vector{Num}, constraint::Interval) → SharedDAG - return ConstraintDAG(var_nodes, all_nodes, root, constraint) +Convenience: build a SharedDAG with a single expression and constraint. +""" +function build_dag(expr::Num, vars::Vector{Num}, constraint::Interval) + dag = SharedDAG(vars) + add_expression!(dag, expr, constraint) + return dag end """ -Recursively build DAG nodes from a symbolic expression. -Returns the node corresponding to `expr`. +Recursively build DAG nodes from a symbolic expression, reusing existing +nodes in the SharedDAG via the node_cache (CSE). """ -function _build_node!(expr, var_map, cache, all_nodes) - # Check cache first (CSE) +function _build_node!(expr, dag::SharedDAG) h = objectid(expr) - haskey(cache, h) && return cache[h] - - node = if haskey(var_map, expr) - # It's a known variable - vn = var_map[expr] - if vn ∉ all_nodes - push!(all_nodes, vn) + haskey(dag.node_cache, h) && return dag.node_cache[h] + + node = if Symbolics.issym(expr) + name = Symbolics.nameof(expr) + if haskey(dag.var_map, name) + dag.var_map[name] + else + error("Unknown variable in expression: $name") end - vn - - elseif Symbolics.issym(expr) - # Unknown symbol — shouldn't happen if vars are correctly specified - error("Unknown variable in expression: $(Symbolics.nameof(expr))") elseif SU.iscall(expr) - _build_operation!(expr, var_map, cache, all_nodes) + _build_operation!(expr, dag) else # Literal number val = Symbolics.value(expr) cn = ConstantNode(val) - push!(all_nodes, cn) + push!(dag.nodes, cn) cn end - cache[h] = node + dag.node_cache[h] = node return node end """ -Build an OperationNode from a symbolic function call. -Handles n-ary `+` and `*` by decomposing to left-associated binary chains. +Build an OperationNode, decomposing n-ary `+` and `*` into binary chains. """ -function _build_operation!(expr, var_map, cache, all_nodes) +function _build_operation!(expr, dag::SharedDAG) op = Symbolics.operation(expr) args = Symbolics.arguments(expr) if (op === (+) || op === (*)) && length(args) > 2 - # Decompose n-ary to binary: (((a + b) + c) + d) - node = _build_node!(args[1], var_map, cache, all_nodes) + node = _build_node!(args[1], dag) for i in 2:length(args) - right = _build_node!(args[i], var_map, cache, all_nodes) + right = _build_node!(args[i], dag) parent = OperationNode(op, AbstractNode[node, right]) - push!(all_nodes, parent) + push!(dag.nodes, parent) node = parent end return node end - # Standard case: unary or binary operation - child_nodes = AbstractNode[ - _build_node!(a, var_map, cache, all_nodes) for a in args - ] - + child_nodes = AbstractNode[_build_node!(a, dag) for a in args] node = OperationNode(op, child_nodes) - push!(all_nodes, node) + push!(dag.nodes, node) return node end diff --git a/src/dag/contractor.jl b/src/dag/contractor.jl index 6ab9a40..96347ab 100644 --- a/src/dag/contractor.jl +++ b/src/dag/contractor.jl @@ -1,30 +1,46 @@ """ DAGContractor <: AbstractContractor -A contractor based on explicit DAG propagation (HC4Revise). -Alternative to the code-generation approach in `Contractor`. - -Stores the expression and variables so it can build a DAG with any constraint -at call time. By default contracts onto `f(x) = 0` (the boundary), matching -the behavior of the code-generation `Contractor`. +A contractor based on persistent DAG propagation (HC4Revise). +The DAG is built once and reused across calls. When multiple contractors +share the same `SharedDAG`, propagating one tightens bounds that benefit +all others via shared variable nodes and subexpressions. + +Constructor forms: +- `DAGContractor(expr, vars)` — builds its own DAG +- `DAGContractor(dag, expr, vars)` — uses an existing SharedDAG (adds the expression to it) """ struct DAGContractor{V} <: AbstractContractor vars::V + dag::SharedDAG ex::Any end -DAGContractor(ex::Num, vars::Vector{Num}) = DAGContractor(vars, ex) +function DAGContractor(ex::Num, vars::Vector{Num}) + dag = SharedDAG(vars) + add_expression!(dag, ex, interval(0.0)) + return DAGContractor(vars, dag, ex) +end + +function DAGContractor(dag::SharedDAG, ex::Num, vars::Vector{Num}) + add_expression!(dag, ex, interval(0.0)) + return DAGContractor(vars, dag, ex) +end function (C::DAGContractor)(X::IntervalBox, constraint::Interval=interval(0.0)) - dag = build_dag(C.ex, C.vars, constraint) - return propagate!(dag, X) + return propagate!(C.dag, X, constraint) end """ DAGSeparator <: AbstractSeparator -A separator that uses DAG-based propagation. -Drop-in replacement for `Separator` — same `(boundary, inner, outer)` output. +A separator that uses persistent DAG-based propagation. +Drop-in replacement for `Separator`. + +Constructor forms: +- `DAGSeparator(expr, vars, constraint)` — builds its own DAG +- `DAGSeparator(dag, expr, vars, constraint)` — uses an existing SharedDAG +- `DAGSeparator(symbolic_inequality, vars)` — parses the inequality """ struct DAGSeparator{V, E, C, F} <: AbstractSeparator vars::V @@ -40,17 +56,21 @@ function DAGSeparator(expr::Num, vars::Vector{Num}, constraint::Interval) return DAGSeparator(vars, expr, constraint, f, contractor) end +function DAGSeparator(dag::SharedDAG, expr::Num, vars::Vector{Num}, constraint::Interval) + contractor = DAGContractor(dag, expr, vars) + f = make_function(expr, vars) + return DAGSeparator(vars, expr, constraint, f, contractor) +end + function (SS::DAGSeparator)(X) if any(x -> isinf(diam(x)), X) return _separate_infinite_dag(SS, X) end - # Contract onto the boundary (f(x) = 0), same as original Separator boundary = SS.contractor(X) lb = IntervalBox(inf.(X)) ub = IntervalBox(sup.(X)) - inner = boundary outer = boundary @@ -74,7 +94,6 @@ end function _separate_infinite_dag(SS::DAGSeparator, X::IntervalBox) C = SS.contractor a, b = inf(SS.constraint), sup(SS.constraint) - inner = C(X, interval(a, b)) local outer @@ -83,9 +102,7 @@ function _separate_infinite_dag(SS::DAGSeparator, X::IntervalBox) elseif b == Inf outer = C(X, interval(-Inf, a)) else - outer1 = C(X, interval(-Inf, a)) - outer2 = C(X, interval(b, Inf)) - outer = outer1 ⊔ outer2 + outer = C(X, interval(-Inf, a)) ⊔ C(X, interval(b, Inf)) end boundary = inner ⊓ outer diff --git a/src/dag/nodes.jl b/src/dag/nodes.jl index 8f577c4..e399d5f 100644 --- a/src/dag/nodes.jl +++ b/src/dag/nodes.jl @@ -21,8 +21,6 @@ end OperationNode(op, children) = OperationNode(op, children, entireinterval()) -# Convenience accessors - is_variable(n::VariableNode) = true is_variable(::AbstractNode) = false @@ -45,23 +43,52 @@ function Base.show(io::IO, n::OperationNode) end """ - ConstraintDAG +A root-constraint pair: an expression root node together with the interval +constraint that its output must satisfy. +""" +struct ConstraintEntry + root::AbstractNode + constraint::Interval +end -Directed acyclic graph representing an expression with a constraint on its output. +""" + SharedDAG -- `variables`: input variable nodes (ordered to match IntervalBox) -- `nodes`: all nodes in topological order (leaves first, root last) -- `root`: the output node -- `constraint`: interval the root must lie in +A persistent DAG that can hold multiple expressions sharing variable nodes +and common subexpressions. Each expression is registered as a `ConstraintEntry` +(root node + constraint interval). + +Built incrementally via `add_expression!`. The `nodes` vector is kept in +topological order (leaves first), so forward propagation is a single pass +and backward propagation is the reverse. """ -struct ConstraintDAG +mutable struct SharedDAG variables::Vector{VariableNode} - nodes::Vector{AbstractNode} # topological order - root::AbstractNode - constraint::Interval + nodes::Vector{AbstractNode} # topological order + constraints::Vector{ConstraintEntry} + var_map::Dict{Symbol, VariableNode} + node_cache::Dict{UInt, AbstractNode} +end + +function SharedDAG(variables::Vector{Num}) + unwrapped = Symbolics.value.(variables) + var_nodes = VariableNode[] + var_map = Dict{Symbol, VariableNode}() + all_nodes = AbstractNode[] + + for (i, v) in enumerate(unwrapped) + name = Symbolics.nameof(v) + vn = VariableNode(name, i, entireinterval()) + push!(var_nodes, vn) + push!(all_nodes, vn) + var_map[name] = vn + end + + return SharedDAG(var_nodes, all_nodes, ConstraintEntry[], var_map, Dict{UInt, AbstractNode}()) end -function Base.show(io::IO, dag::ConstraintDAG) +function Base.show(io::IO, dag::SharedDAG) vars = join([string(v.name) for v in dag.variables], ", ") - print(io, "ConstraintDAG($(length(dag.nodes)) nodes, vars=($vars), constraint=$(dag.constraint))") + nc = length(dag.constraints) + print(io, "SharedDAG($(length(dag.nodes)) nodes, vars=($vars), $nc constraints)") end diff --git a/src/dag/propagate.jl b/src/dag/propagate.jl index 50511bd..9e66021 100644 --- a/src/dag/propagate.jl +++ b/src/dag/propagate.jl @@ -1,22 +1,18 @@ """ - forward!(dag::ConstraintDAG, X::IntervalBox) + forward!(dag::SharedDAG, X::IntervalBox) -Forward propagation: evaluate intervals from leaves to root in topological order. -Sets `range` on every node. +Forward propagation: set variable ranges from X, then evaluate all operation +nodes in topological order. """ -function forward!(dag::ConstraintDAG, X::IntervalBox) +function forward!(dag::SharedDAG, X::IntervalBox) for v in dag.variables v.range = X[v.index] end - for node in dag.nodes if node isa OperationNode node.range = _eval_forward(node.op, node.children) end - # VariableNode and ConstantNode ranges are already set end - - return dag.root.range end function _eval_forward(op::Function, children::Vector{AbstractNode}) @@ -25,13 +21,10 @@ function _eval_forward(op::Function, children::Vector{AbstractNode}) elseif length(children) == 2 return op(children[1].range, children[2].range) else - # Shouldn't happen after n-ary decomposition, but handle gracefully return op([c.range for c in children]...) end end -# Special case: ^ with integer exponent. Symbolics stores ^(x, 2) where -# the second arg is a ConstantNode. IntervalArithmetic dispatches x^n. function _eval_forward(::typeof(^), children::Vector{AbstractNode}) base = children[1].range exp_node = children[2] @@ -44,28 +37,28 @@ function _eval_forward(::typeof(^), children::Vector{AbstractNode}) return base ^ exp_node.range end - """ - backward!(dag::ConstraintDAG) + backward!(dag::SharedDAG) -Backward propagation: propagate constraint from root to leaves in reverse -topological order using reverse interval contractors. +Backward propagation for all constraints. For each constraint entry, intersect +the root's range with the constraint, then propagate backward through the DAG. -The root's range is first intersected with the constraint. Then each -OperationNode propagates narrowed ranges to its children. +When multiple constraints share nodes, each backward pass further narrows +shared variable and subexpression ranges. """ -function backward!(dag::ConstraintDAG) - dag.root.range = IntervalArithmetic.intersect_interval(dag.root.range, dag.constraint) +function backward!(dag::SharedDAG) + for entry in dag.constraints + entry.root.range = IntervalArithmetic.intersect_interval(entry.root.range, entry.constraint) + end - if isempty_interval(dag.root.range) - # Constraint is infeasible on this box — mark everything empty + all_empty = all(isempty_interval(e.root.range) for e in dag.constraints) + if all_empty for v in dag.variables v.range = emptyinterval() end return end - # Traverse in reverse topological order (root → leaves) for i in length(dag.nodes):-1:1 node = dag.nodes[i] node isa OperationNode || continue @@ -74,50 +67,67 @@ function backward!(dag::ConstraintDAG) end """ -Apply the reverse contractor for a single operation node, narrowing -both `node.range` and each child's range in place. + backward!(dag::SharedDAG, constraint::Interval) + +Backward propagation with a single override constraint applied to the first +(or only) root. Used by DAGContractor when contracting with a specific constraint. """ +function backward!(dag::SharedDAG, constraint::Interval) + entry = dag.constraints[1] + entry.root.range = IntervalArithmetic.intersect_interval(entry.root.range, constraint) + + if isempty_interval(entry.root.range) + for v in dag.variables + v.range = emptyinterval() + end + return + end + + for i in length(dag.nodes):-1:1 + node = dag.nodes[i] + node isa OperationNode || continue + _contract_backward!(node) + end +end + function _contract_backward!(node::OperationNode) - op = node.op - children = node.children z = node.range - + children = node.children if isempty_interval(z) for c in children c.range = emptyinterval() end return end - if length(children) == 1 - _contract_unary!(op, z, children[1]) + _contract_unary!(node.op, z, children[1]) elseif length(children) == 2 - _contract_binary!(op, z, children[1], children[2]) + _contract_binary!(node.op, z, children[1], children[2]) end end # --- Binary reverse contractors --- function _contract_binary!(::typeof(+), z, left, right) - z_new, l_new, r_new = plus_rev(z, left.range, right.range) + _, l_new, r_new = plus_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end function _contract_binary!(::typeof(-), z, left, right) - z_new, l_new, r_new = minus_rev(z, left.range, right.range) + _, l_new, r_new = minus_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end function _contract_binary!(::typeof(*), z, left, right) - z_new, l_new, r_new = mul_rev(z, left.range, right.range) + _, l_new, r_new = mul_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end function _contract_binary!(::typeof(/), z, left, right) - z_new, l_new, r_new = div_rev(z, left.range, right.range) + _, l_new, r_new = div_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end @@ -126,26 +136,24 @@ function _contract_binary!(::typeof(^), z, base_node, exp_node) if exp_node isa ConstantNode n = exp_node.value if n == round(Int, n) - # power_rev(z, x, n::Integer) → (z', x', n) - z_new, x_new, _ = power_rev(z, base_node.range, round(Int, n)) + _, x_new, _ = power_rev(z, base_node.range, round(Int, n)) base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new) return end end - # General case: x^y - z_new, x_new, y_new = power_rev(z, base_node.range, exp_node.range) + _, x_new, y_new = power_rev(z, base_node.range, exp_node.range) base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new) exp_node.range = IntervalArithmetic.intersect_interval(exp_node.range, y_new) end function _contract_binary!(::typeof(min), z, left, right) - z_new, l_new, r_new = min_rev(z, left.range, right.range) + _, l_new, r_new = min_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end function _contract_binary!(::typeof(max), z, left, right) - z_new, l_new, r_new = max_rev(z, left.range, right.range) + _, l_new, r_new = max_rev(z, left.range, right.range) left.range = IntervalArithmetic.intersect_interval(left.range, l_new) right.range = IntervalArithmetic.intersect_interval(right.range, r_new) end @@ -153,43 +161,27 @@ end # --- Unary reverse contractors --- for (jl_op, rev_fn) in [ - (:sin, :sin_rev), - (:cos, :cos_rev), - (:tan, :tan_rev), - (:asin, :asin_rev), - (:acos, :acos_rev), - (:atan, :atan_rev), - (:sinh, :sinh_rev), - (:cosh, :cosh_rev), - (:tanh, :tanh_rev), - (:asinh, :asinh_rev), - (:acosh, :acosh_rev), - (:atanh, :atanh_rev), - (:exp, :exp_rev), - (:exp2, :exp2_rev), - (:exp10, :exp10_rev), + (:sin, :sin_rev), (:cos, :cos_rev), (:tan, :tan_rev), + (:asin, :asin_rev), (:acos, :acos_rev), (:atan, :atan_rev), + (:sinh, :sinh_rev), (:cosh, :cosh_rev), (:tanh, :tanh_rev), + (:asinh, :asinh_rev), (:acosh, :acosh_rev), (:atanh, :atanh_rev), + (:exp, :exp_rev), (:exp2, :exp2_rev), (:exp10, :exp10_rev), (:expm1, :expm1_rev), - (:log, :log_rev), - (:log2, :log2_rev), - (:log10, :log10_rev), + (:log, :log_rev), (:log2, :log2_rev), (:log10, :log10_rev), (:log1p, :log1p_rev), - (:sqrt, :sqrt_rev), - (:abs, :abs_rev), - (:sign, :sign_rev), + (:sqrt, :sqrt_rev), (:abs, :abs_rev), (:sign, :sign_rev), ] @eval function _contract_unary!(::typeof($jl_op), z, child) - z_new, x_new = $rev_fn(z, child.range) + _, x_new = $rev_fn(z, child.range) child.range = IntervalArithmetic.intersect_interval(child.range, x_new) end end -# Unary minus: -(x) function _contract_unary!(::typeof(-), z, child) - z_new, x_new = minus_rev(z, child.range) + _, x_new = minus_rev(z, child.range) child.range = IntervalArithmetic.intersect_interval(child.range, x_new) end -# Fallback for unknown operations — no contraction, just leave ranges as-is function _contract_unary!(op, z, child) @warn "No reverse contractor for unary $op — skipping" maxlog=1 end @@ -198,42 +190,42 @@ function _contract_binary!(op, z, left, right) @warn "No reverse contractor for binary $op — skipping" maxlog=1 end - """ - propagate!(dag::ConstraintDAG, X::IntervalBox; max_iter=10, tol=1e-10) → IntervalBox - -Iterate forward-backward propagation on the DAG until convergence or `max_iter` -iterations. Returns the narrowed IntervalBox. + propagate!(dag::SharedDAG, X::IntervalBox; max_iter=10, tol=1e-10) → IntervalBox -This implements the HC4Revise algorithm: each iteration does one forward pass -(natural interval evaluation) followed by one backward pass (constraint narrowing). -For expressions where a variable appears multiple times, a single pass cannot -fully propagate — iteration is needed because each pass tightens the shared node, -which then feeds tighter bounds into the next forward evaluation. +Iterate forward-backward propagation on the DAG until convergence. +All constraints in the DAG are propagated each iteration, so narrowing +from one constraint benefits the others via shared variable nodes. """ -function propagate!(dag::ConstraintDAG, X::IntervalBox; max_iter::Int=10, tol::Float64=1e-10) +function propagate!(dag::SharedDAG, X::IntervalBox; max_iter::Int=10, tol::Float64=1e-10) for _ in 1:max_iter old_ranges = IntervalBox(Tuple(v.range for v in dag.variables)) - forward!(dag, X) backward!(dag) - - # Build new X from narrowed variable ranges X = IntervalBox(Tuple(v.range for v in dag.variables)) + isempty(X) && return X + max_change = maximum(diam(old) - diam(new) for (old, new) in zip(old_ranges, X)) + max_change < tol && break + end + return X +end - if isempty(X) - return X - end +""" + propagate!(dag::SharedDAG, X::IntervalBox, constraint::Interval; ...) → IntervalBox - # Check convergence: did the variable ranges change significantly? - max_change = maximum( - diam(old) - diam(new) - for (old, new) in zip(old_ranges, X) - ) - if max_change < tol - break - end +Propagate with a single override constraint (ignoring stored constraints). +Used by DAGContractor for boundary/complement contraction. +""" +function propagate!(dag::SharedDAG, X::IntervalBox, constraint::Interval; + max_iter::Int=10, tol::Float64=1e-10) + for _ in 1:max_iter + old_ranges = IntervalBox(Tuple(v.range for v in dag.variables)) + forward!(dag, X) + backward!(dag, constraint) + X = IntervalBox(Tuple(v.range for v in dag.variables)) + isempty(X) && return X + max_change = maximum(diam(old) - diam(new) for (old, new) in zip(old_ranges, X)) + max_change < tol && break end - return X end diff --git a/test/test_dag.jl b/test/test_dag.jl index cbf62ee..2ca8c48 100644 --- a/test/test_dag.jl +++ b/test/test_dag.jl @@ -1,67 +1,93 @@ using IntervalArithmetic, IntervalArithmetic.Symbols using IntervalConstraintProgramming -using IntervalConstraintProgramming: build_dag, forward!, backward!, propagate!, ConstraintDAG +using IntervalConstraintProgramming: build_dag, forward!, backward!, propagate!, + SharedDAG, add_expression!, ConstraintEntry using Symbolics using IntervalBoxes using Test -const IntervalType{T} = Union{Interval{T}, BareInterval{T}} -eq(a::IntervalType, b::IntervalType) = isequal_interval(bareinterval(a), bareinterval(b)) +const IType{T} = Union{Interval{T}, BareInterval{T}} +eq(a::IType, b::IType) = isequal_interval(bareinterval(a), bareinterval(b)) eq(a::IntervalBox, b::IntervalBox) = all(eq.(a, b)) eq(a::Vector, b::Vector) = length(a) == length(b) && all(eq.(a, b)) -@testset "DAG construction" begin +@testset "SharedDAG construction" begin vars = @variables x, y - dag = build_dag(x + y, vars, interval(0.0)) + dag = SharedDAG(vars) @test length(dag.variables) == 2 + @test length(dag.nodes) == 2 # just variable nodes @test dag.variables[1].name == :x @test dag.variables[2].name == :y - dag2 = build_dag(x^2 + y^2, vars, interval(-Inf, 0.0)) - # Should have: x, y, 2 (const), x^2, y^2, x^2+y^2 = 6 nodes min - @test length(dag2.nodes) >= 6 + add_expression!(dag, x^2 + y^2, interval(-Inf, 0.0)) + @test length(dag.constraints) == 1 + @test length(dag.nodes) > 2 +end + +@testset "SharedDAG node reuse (CSE)" begin + vars = @variables x, y + + dag = SharedDAG(vars) + add_expression!(dag, x^2 + y^2 - 1, interval(-Inf, 0.0)) + add_expression!(dag, x^2 - y, interval(-Inf, 0.0)) + + pow_nodes = filter(n -> n isa IntervalConstraintProgramming.OperationNode && n.op === (^), dag.nodes) + @test length(pow_nodes) == 2 # x^2 and y^2, not a duplicate x^2 + + @test length(dag.constraints) == 2 +end + +@testset "SharedDAG n-ary decomposition" begin + vars = @variables x, y - # N-ary: x + y + 1 should decompose to binary - dag3 = build_dag(x + y + 1, vars, interval(0.0)) - ops = filter(n -> n isa IntervalConstraintProgramming.OperationNode, dag3.nodes) - # Each + node should have exactly 2 children + dag = build_dag(x + y + 1, vars, interval(0.0)) + ops = filter(n -> n isa IntervalConstraintProgramming.OperationNode, dag.nodes) for op in ops @test length(op.children) <= 2 end end -@testset "DAG forward propagation" begin +@testset "Forward propagation" begin vars = @variables x, y dag = build_dag(x^2 + y^2, vars, interval(-Inf, 1.0)) X = IntervalBox(interval(0.0, 1.0), interval(0.0, 1.0)) - result = forward!(dag, X) - # x^2 ∈ [0, 1], y^2 ∈ [0, 1], sum ∈ [0, 2] but with the -1 from analyse... - # Actually analyse(x^2+y^2 <= 1) gives expr = -1 + x^2 + y^2, constraint = (-∞, 0] - # Let's test the raw expression directly - @test inf(result) <= 2.0 - @test sup(result) >= 0.0 + forward!(dag, X) + root_range = dag.constraints[1].root.range + @test inf(root_range) <= 0.0 + 1e-10 + @test sup(root_range) >= 2.0 - 1e-10 end -@testset "DAG propagation contracts correctly" begin +@testset "Single-expression propagation" begin vars = @variables x, y - # x^2 + y^2 - 1 ≤ 0 (the unit disk) ex, constraint = IntervalConstraintProgramming.analyse(x^2 + y^2 <= 1) dag = build_dag(ex, vars, constraint) X = IntervalBox(interval(-10.0, 10.0), 2) result = propagate!(dag, X) - # Should narrow to [-1, 1]² @test all(inf.(result) .>= -1.0 - 1e-10) @test all(sup.(result) .<= 1.0 + 1e-10) +end + +@testset "Multi-constraint propagation" begin + vars = @variables x, y - # Box already inside: should not widen - X2 = IntervalBox(interval(0.0, 0.3), 2) - dag2 = build_dag(ex, vars, constraint) - result2 = propagate!(dag2, X2) - @test all(inf.(result2) .>= 0.0 - 1e-10) - @test all(sup.(result2) .<= 0.3 + 1e-10) + dag = SharedDAG(vars) + add_expression!(dag, x^2 + y^2 - 1, interval(-Inf, 0.0)) + add_expression!(dag, x - y, interval(0.0, Inf)) + + X = IntervalBox(interval(-10, 10), 2) + result = propagate!(dag, X) + @test all(inf.(result) .>= -1.0 - 1e-10) + @test all(sup.(result) .<= 1.0 + 1e-10) + + # y ≥ x² narrows y ≥ 0 + dag2 = SharedDAG(vars) + add_expression!(dag2, x^2 + y^2 - 1, interval(-Inf, 0.0)) + add_expression!(dag2, x^2 - y, interval(-Inf, 0.0)) + result2 = propagate!(dag2, IntervalBox(interval(-10, 10), 2)) + @test inf(result2[2]) >= -1e-10 # y ≥ 0 end @testset "DAGContractor" begin @@ -72,19 +98,27 @@ end result = C(X, -Inf..1) @test eq(result, IntervalBox(-1..1, 2)) - X2 = IntervalBox(interval(0.5, 1.5), 3) vars3 = @variables x, y, z + X3 = IntervalBox(interval(0.5, 1.5), 3) C1 = DAGContractor(x + y, vars3) - @test eq(C1(X2, -Inf..1), IntervalBox(interval(0.5, 0.5), interval(0.5, 0.5), interval(0.5, 1.5))) + @test eq(C1(X3, -Inf..1), IntervalBox(interval(0.5, 0.5), interval(0.5, 0.5), interval(0.5, 1.5))) +end - C2 = DAGContractor(y + z, vars3) - @test eq(C2(X2, -Inf..1), IntervalBox(interval(0.5, 1.5), interval(0.5, 0.5), interval(0.5, 0.5))) +@testset "DAGContractor with shared DAG" begin + vars = @variables x, y + + dag = SharedDAG(vars) + C1 = DAGContractor(dag, x^2 + y^2, vars) + C2 = DAGContractor(dag, x - y, vars) + + @test C1.dag === C2.dag + @test length(dag.constraints) == 2 end @testset "DAGSeparator" begin vars = @variables x, y - II = -100..100 + II = interval(-100, 100) X = IntervalBox(II, 2) S = DAGSeparator(x^2 + y^2 <= 1, vars) @@ -98,19 +132,17 @@ end @test eq(outer2, IntervalBox(-Inf..Inf, 2)) end -@testset "DAG pave" begin +@testset "DAGSeparator with shared DAG" begin vars = @variables x, y - S = DAGSeparator(x^2 + y^2 <= 1, vars) - X = IntervalBox(-3..3, 2) - inner, boundary = pave(X, S, 0.5) - @test length(inner) > 0 - @test length(boundary) > 0 - - inner_area = sum(prod(diam.(b)) for b in inner) - total_area = inner_area + sum(prod(diam.(b)) for b in boundary) - @test inner_area > 2.0 # should be close to π - @test total_area < 5.0 # should bracket π from above + dag = SharedDAG(vars) + ex1, c1 = IntervalConstraintProgramming.analyse(x^2 + y^2 <= 1) + ex2, c2 = IntervalConstraintProgramming.analyse(x >= 0) + S1 = DAGSeparator(dag, ex1, vars, c1) + S2 = DAGSeparator(dag, ex2, vars, c2) + + @test S1.contractor.dag === S2.contractor.dag + @test length(dag.constraints) == 2 end @testset "DAG pave matches original" begin @@ -119,7 +151,7 @@ end S_orig = Separator(x^2 + y^2 <= 1, vars) S_dag = DAGSeparator(x^2 + y^2 <= 1, vars) - X = IntervalBox(-3..3, 2) + X = IntervalBox(interval(-3, 3), 2) inner_orig, boundary_orig = pave(X, S_orig, 0.5) inner_dag, boundary_dag = pave(X, S_dag, 0.5) @@ -134,7 +166,7 @@ end S2 = DAGSeparator(y > 0, vars) S = S1 ⊓ S2 - X = IntervalBox(-3..3, 2) + X = IntervalBox(interval(-3, 3), 2) boundary, inner, outer = S(X) @test eq(inner, IntervalBox(0..3, 2)) end @@ -143,7 +175,7 @@ end vars = @variables x, y, z S = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, vars) - X = IntervalBox(-10..10, 3) + X = IntervalBox(interval(-10, 10), 3) inner, boundary = pave(X, S, 1.0) @test inner[1] isa IntervalBox{3, Float64, Interval{Float64}} end