Compiled DAG propagation: ~3x speedup over persistent DAG

.github/workflows/CI.yml

@@ -12,6 +12,7 @@ jobs:
       matrix:
         version:
           - '1.10'
+          - '1.12'
           - 'nightly'
         os:
           - ubuntu-latest
@@ -24,12 +25,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 +43,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

.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 }}

CLAUDE.md

@@ -0,0 +1,88 @@
+# 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`.
+
+## DAG-based constraint propagation (2026-04-14)
+
+### Overview
+
+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/`.
+
+**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.
+
+**Files:** `src/dag/{nodes,build,propagate,contractor}.jl`, `test/test_dag.jl`, `benchmark/bench_dag_vs_codegen.jl`
+
+### Architecture
+
+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.
+
+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.
+
+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.
+
+### Performance (benchmarks, persistent shared DAG)
+
+| Problem | ϵ | Codegen | DAG (persistent) | Ratio |
+|---------|---|---------|-----------------|-------|
+| Unit disk 2D | single call | 1.2 μs | 6.8 μs | 6x |
+| Unit disk 2D | 0.1 | 294 μs | 1.0 ms | 3.5x |
+| Unit disk 2D | 0.01 | 2.5 ms | 8.8 ms | 3.5x |
+| Annulus 2D | 0.1 | 1.6 ms | 4.4 ms | 2.8x |
+| 3D torus | 1.0 | 12 ms | 33 ms | 2.8x |
+
+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]`.
+
+### Compiled propagation schedule
+
+The `CompiledDAG` (`src/dag/compile.jl`) converts the SharedDAG into a flat instruction array operating on a pre-allocated workspace of intervals. This combines three optimizations:
+
+1. **Symbol-keyed operations** — `op::Symbol` (`:add`, `:mul`, etc.) instead of `op::Function`, enabling efficient `===` comparison without dynamic dispatch
+2. **Pre-allocated workspace** — flat `Vector{Interval}` indexed by slot number, no per-node allocation during traversal
+3. **Flat instruction array** — `Vector{Instruction}` in topological order; forward iterates forward, backward iterates in reverse. Each `Instruction` stores `(op, out_slot, arg1_slot, arg2_slot)`.
+
+The compiled version is ~3x faster than the uncompiled persistent DAG and ~3x slower than the code-generation approach.

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"
@@ -11,13 +11,13 @@ 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 = "0.22.12 - 0.23, 1"
+IntervalBoxes = "0.2 - 0.3"
+IntervalContractors = "0.5 - 0.6"
+ReversePropagation = "0.3 - 0.4"
 StaticArrays = "1"
-Symbolics = "5, 6"
-julia = "1"
+Symbolics = "5 - 7"
+julia = "1.10"
 
 [extras]
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

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

benchmark/bench_dag_vs_codegen.jl

@@ -0,0 +1,107 @@
+using IntervalConstraintProgramming
+using IntervalConstraintProgramming: SharedDAG, add_expression!, propagate!
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalBoxes
+using Symbolics
+using BenchmarkTools
+
+@variables x, y, z
+
+println("=" ^ 60)
+println("Benchmark: DAG (persistent) vs code-generation")
+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(interval(-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(interval(-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(interval(-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")
+
+# --- 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))