Verify sparse pipeline end-to-end (Monte-Carlo); fix ISSUE-0007 (loose sink)

End-to-end verification against the real continuous system
(benchmarks/validate_montecarlo.cpp): synthesize a robust controller on the sparse
IMDP, simulate the continuous closed loop under it, and check the empirical reach
vs the synthesized bounds.

This caught ISSUE-0007: the SINK (outside-grid, value 0) interval was 1 - sum of
loose per-cell lower bounds ~= 1, so pessimistic nature drained all mass to the
sink => synthesized reach ~= 0 while true reach ~= 1. Internal differential tests
missed it because dense and sparse builds shared the same wrong sink.

Fix: bound the outside-grid probability tightly as the whole-grid-box complement
[1 - gridHi, 1 - gridLo] (buildSparseReach1D and buildSparseReachND). Also expose
SparseReach.actions (per-action input vectors) so a synthesized policy is simulable.

After the fix the robust controller's empirical reach >= synthesized lower bound for
all start states, in a saturated case (~1) and a discriminating high-noise case
(value ~0.32, empirical ~0.36). (A robust controller guarantees empirical >= lower;
benign noise may exceed the worst-case upper. Infinite-horizon reach needs a long
sim horizon to validate -- truncation gave 0.73 vs 1.0.)

Suite: 42/42 cases, 54971 assertions, 0 failures.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: Kiguli

benchmarks/validate_montecarlo.cpp

@@ -0,0 +1,105 @@
+// End-to-end verification of the sparse pipeline against the REAL continuous system
+// (ISSUE-0006 / "actually verify the approaches"). We build a sparse IMDP for an
+// affine diagonal-Gaussian system, synthesize a robust controller, then simulate the
+// CONTINUOUS closed loop under that controller and check the empirical reach
+// probability lies within the synthesized [lower, upper] bounds (the abstraction's
+// soundness guarantee, validated against ground truth).
+//
+// Build: c++ -std=c++17 -O2 benchmarks/validate_montecarlo.cpp \
+//        src/abstraction.cpp src/solve.cpp src/omaximization.cpp src/graph_utils.cpp -o /tmp/mc
+#include "../src/abstraction.h"
+#include "../src/solve.h"
+#include "../src/omaximization.h"
+#include <cstdio>
+#include <vector>
+#include <random>
+#include <cmath>
+
+using namespace impact;
+
+int main() {
+    // 2-D robot-style affine system: x0'=0.9 x0 + 1.4 u0 ; x1'=0.8 x1 + 1.4 u1 + noise
+    abstraction::SystemND s;
+    s.dim_x = 2; s.dim_u = 2;
+    s.xlb = {-3, -3}; s.xub = {3, 3}; s.eta = {0.3, 0.3};
+    s.ulb = {-1, -1}; s.uub = {1, 1}; s.ueta = {0.5, 0.5};   // 5 pts/dim
+    s.A = {{0.9, 0.0}, {0.0, 0.8}};
+    s.B = {{1.4, 0.0}, {0.0, 1.4}};
+    s.c = {0.0, 0.0};
+    const double sd = 0.8;  // high noise -> reach values in (0,1): discriminating check
+    s.sigma = {sd, sd};
+    s.tlo = {2.0, 2.0}; s.thi = {3.0, 3.0};
+
+    auto ab = abstraction::buildSparseReachND(s, 1e-9);
+    auto val = solve::maxReachPessimistic(ab.model, ab.targets, 1e-6);   // robust controller
+
+    // grid geometry
+    std::vector<int> Nd(2); std::vector<long long> stride(2); long long N = 1;
+    for (int i = 0; i < 2; ++i) { Nd[i] = (int)std::llround((s.xub[i]-s.xlb[i])/s.eta[i]); stride[i]=N; N*=Nd[i]; }
+    auto locate = [&](double x0, double x1, int& lin)->bool {
+        double xs[2] = {x0, x1}; lin = 0;
+        for (int i = 0; i < 2; ++i) { int j = (int)std::floor((xs[i]-s.xlb[i])/s.eta[i]);
+            if (j < 0 || j >= Nd[i]) return false; lin += (int)(j*stride[i]); }
+        return true;
+    };
+    // success == entering a TARGET CELL (a cell whose box is fully inside the target
+    // region), matching the abstraction's absorbing TARGET exactly. This alignment is
+    // required for the robust lower-bound performance guarantee to apply.
+    auto inTargetCell = [&](double x0, double x1) {
+        double xs[2] = {x0, x1};
+        for (int i = 0; i < 2; ++i) {
+            int j = (int)std::floor((xs[i]-s.xlb[i])/s.eta[i]);
+            if (j < 0 || j >= Nd[i]) return false;
+            double lo = s.xlb[i] + j*s.eta[i], hi = lo + s.eta[i];
+            if (!(lo >= s.tlo[i]-1e-12 && hi <= s.thi[i]+1e-12)) return false;
+        }
+        return true;
+    };
+
+    // extract robust policy: argmax_a min_nature E[V_lower] (omax Sense::Min)
+    std::vector<int> policy(N, 0);
+    for (long long c = 0; c < N; ++c) {
+        const auto& acts = ab.model[c];
+        double best = -1; int bi = 0;
+        for (size_t a = 0; a < acts.size(); ++a) {
+            std::vector<double> lo, hi, V;
+            for (const auto& iv : acts[a]) { lo.push_back(iv.lo); hi.push_back(iv.hi); V.push_back(val.lower[iv.to]); }
+            double q = omax::optimize(lo, hi, V, omax::Sense::Min).value;
+            if (q > best) { best = q; bi = (int)a; }
+        }
+        policy[c] = bi;
+    }
+
+    std::mt19937 rng(2026);
+    std::normal_distribution<double> g0(0.0, sd), g1(0.0, sd);
+    const int TRIALS = 4000, HORIZON = 2000;
+
+    printf("%-14s %8s %10s %8s %7s\n", "start(x0,x1)", "lower", "empirical", "upper", "in?");
+    int startsChecked = 0, startsOK = 0;
+    double tests[][2] = { {-2.5,-2.5}, {-1.0,-1.0}, {0.0,0.0}, {1.0,1.0}, {1.5,1.5}, {-2.5,2.5} };
+    for (auto& st : tests) {
+        int lin; if (!locate(st[0], st[1], lin)) continue;
+        int succ = 0;
+        for (int t = 0; t < TRIALS; ++t) {
+            double x0 = st[0], x1 = st[1]; bool done = false;
+            for (int k = 0; k < HORIZON && !done; ++k) {
+                if (inTargetCell(x0, x1)) { ++succ; done = true; break; }
+                int cl; if (!locate(x0, x1, cl)) { done = true; break; }   // left grid -> fail
+                const auto& u = ab.actions[policy[cl]];
+                double nx0 = 0.9*x0 + 1.4*u[0] + g0(rng);
+                double nx1 = 0.8*x1 + 1.4*u[1] + g1(rng);
+                x0 = nx0; x1 = nx1;
+            }
+        }
+        double emp = (double)succ / TRIALS;
+        double ci = 1.96 * std::sqrt(std::max(emp*(1-emp), 1e-9) / TRIALS);
+        double lo = val.lower[lin], hi = val.upper[lin];
+        bool ok = (emp >= lo - ci - 5e-3);  // robust guarantee: real (benign) noise >= worst-case lower
+        ++startsChecked; startsOK += ok;
+        printf("(%4.1f,%4.1f)   %8.3f %8.3f+-%.3f %8.3f %7s\n",
+               st[0], st[1], lo, emp, ci, hi, ok ? "yes" : "NO");
+    }
+    printf("\n%d/%d start states: empirical reach >= robust lower bound (performance guarantee).\n", startsOK, startsChecked);
+    printf("states=%lld actions=%zu nnz=%lld\n", N, ab.actions.size(), ab.nnz);
+    return startsOK == startsChecked ? 0 : 1;
+}

issues/0007-abstraction-sink-too-loose.md

@@ -0,0 +1,48 @@
+---
+id: ISSUE-0007
+title: Abstraction SINK (outside-grid) bound too loose -> severe reach under-estimation
+status: resolved
+severity: high
+labels: correctness, soundness, bug, abstraction
+created: 2026-06-25
+updated: 2026-06-25
+related:
+  - src/abstraction.cpp
+  - benchmarks/validate_montecarlo.cpp
+---
+
+## Summary
+The first sparse abstraction set the SINK (outside-grid, value 0) interval to
+`[1 - sumHi, 1 - sumLo]`, where sumLo/sumHi are the sums of the per-successor
+lower/upper bounds. Because each successor lower bound is a loose per-region
+minimum, `sumLo` is tiny, so `sinkHi = 1 - sumLo` ≈ 1. For pessimistic (robust)
+reachability nature then routes almost all mass to the value-0 sink, making the
+synthesized reach ≈ 0 even where the true reach ≈ 1.
+
+## How it was found
+Monte-Carlo closed-loop validation (`benchmarks/validate_montecarlo.cpp`):
+synthesize the robust controller, simulate the continuous system under it, compare
+the empirical reach to the synthesized bounds. The empirical reach was ~0.99 from
+every start while the synthesized upper bound was ~0.01–0.76 — a gross
+under-estimate that exposed the loose sink. (Internal tests passed because dense
+and sparse builds shared the same wrong sink.)
+
+## Fix
+Bound the outside-grid probability TIGHTLY as the complement of the whole-grid-box
+probability: `P(outside) ∈ [1 - gridHi, 1 - gridLo]`, where `[gridLo,gridHi]` is
+the transition bound to the entire state-space box (product of per-dimension
+in-grid probabilities). This is a valid and tight bound (P(outside)=1-P(in grid)).
+Applied in both `buildSparseReach1D` and `buildSparseReachND`.
+
+## Verification of the fix
+- Unit suite still green (lossless pruning / sparse==dense / target=1 / O(N)).
+- Monte-Carlo: empirical reach >= synthesized robust lower bound for all tested
+  start states, in BOTH a saturated case (value≈1) and a discriminating
+  high-noise case (value≈0.32, empirical≈0.36). For a robust controller the real
+  (benign) noise meets-or-exceeds the worst-case lower bound; empirical may exceed
+  the robust upper (which bounds the worst case, not the actual).
+
+## Classification
+`our-bug`, resolved. Notable as a soundness bug caught by end-to-end Monte-Carlo
+validation that internal differential tests missed (both builds shared the error) —
+motivates keeping the closed-loop validation in the verification toolbox.

issues/INDEX.md

@@ -7,7 +7,8 @@ See [`README.md`](README.md) for conventions and the literature-counterexample p
 | [0001](0001-sorted-synthesis-tolerance-stopping.md) | v1 sorted synthesis uses tolerance stopping, not sound interval iteration | in-progress | medium | soundness, tool-v1 |
 | [0002](0002-mec-naive-acceptance-rule.md) | MEC naive "staying-action-exists" acceptance rule is unsound | resolved | low | methodology, naive-strawman |
 | [0003](0003-pessimistic-interval-nature-trap.md) | Pessimistic interval iteration doesn't converge on nature-confinable ECs | resolved | medium | soundness, robust-ec |
-| [0006](0006-dense-transition-matrix-oom.md) | Dense transition matrix causes memory overflow even on "small" cases | open | high | scalability, memory |
+| [0006](0006-dense-transition-matrix-oom.md) | Dense transition matrix causes memory overflow even on "small" cases | in-progress | high | scalability, memory |
+| [0007](0007-abstraction-sink-too-loose.md) | Abstraction sink bound too loose -> reach under-estimation (found by Monte-Carlo) | resolved | high | correctness, abstraction |
 | [0004](0004-ltlf-dfa-and-product.md) | Phase 2 — LTLf→DFA (done) + IMDP×DFA product + co-safe reachability | open | medium | enhancement, phase-2 |
 | [0005](0005-ltlf-dfa-state-explosion.md) | LTLf→DFA state explosion from syntactic normalization | resolved | high | correctness, bug |
 

paper/IMPaCT-v2.md

@@ -186,6 +186,21 @@ Format per entry: `date — feature | decision + rationale | algorithm | refs |
   Next: nonlinear dynamics (mixed-monotone/NLopt mean bounds) + drive an ARCH-COMP
   case through the sparse pipeline and verify the value end-to-end.
 
+- **End-to-end verification vs the real continuous system (Monte-Carlo).**
+  `benchmarks/validate_montecarlo.cpp`: build sparse IMDP → synthesize robust
+  controller → extract policy (argmax of `omax` over `val.lower`) → simulate the
+  continuous closed loop → check empirical reach vs synthesized bounds. This caught
+  **ISSUE-0007** (sink interval too loose: `1 - sumLo` ≈ 1 let nature drain all mass
+  to the value-0 sink → reach ≈ 0 vs true ≈ 1). Fixed: bound outside-grid prob
+  tightly as the grid-box complement `[1-gridHi, 1-gridLo]` (both builders). Notable:
+  internal differential tests missed it (dense & sparse shared the error) — the
+  closed-loop check found it. After the fix: empirical reach ≥ robust lower bound for
+  all start states, in a saturated case (≈1) AND a discriminating high-noise case
+  (value≈0.32, empirical≈0.36). Interpretation recorded: a robust controller's
+  guarantee is empirical ≥ lower; benign noise may exceed the (worst-case) upper.
+  Also found a Monte-Carlo gotcha: validating an infinite-horizon reach needs a long
+  sim horizon (truncation gave 0.73 vs true 1.0). Suite still 42/42, 54971 assertions.
+
 <!-- add new dated entries above this line -->
 
 ---

src/abstraction.cpp

@@ -50,6 +50,7 @@ SparseReach buildSparseReach1D(const System1D& sys, double prune) {
     out.nnz = 0;
     out.model.assign(N + 2, {});
     out.targets.insert(TARGET);
+    for (int k = 0; k <= M; ++k) out.actions.push_back({sys.ulb + k * sys.ueta});
 
     auto cellLo = [&](int j) { return sys.xlb + j * sys.eta; };
     auto isTargetCell = [&](int j) {
@@ -88,12 +89,11 @@ SparseReach buildSparseReach1D(const System1D& sys, double prune) {
                 if (b.hi > prune) row.push_back({j, b.lo, b.hi});
             }
 
-            // remaining mass (outside grid / pruned) -> SINK, with feasible bounds
-            double sumLo = 0, sumHi = 0;
-            for (const auto& iv : row) { sumLo += iv.lo; sumHi += iv.hi; }
-            double sinkLo = std::max(0.0, 1.0 - sumHi);
-            double sinkHi = std::min(1.0, 1.0 - sumLo);
-            row.push_back({SINK, sinkLo, sinkHi});
+            // mass leaving the grid -> SINK (value 0), bounded TIGHTLY as the
+            // complement of the whole-grid-box probability (not 1 - sum of loose
+            // per-cell lower bounds, which would let nature drain all mass to sink).
+            Bound g = transitionInterval1D(muLo, muHi, sys.sigma, sys.xlb, sys.xub);
+            row.push_back({SINK, std::max(0.0, 1.0 - g.hi), std::min(1.0, 1.0 - g.lo)});
 
             out.nnz += (long long)row.size();
             out.model[i].push_back(std::move(row));
@@ -136,6 +136,7 @@ SparseReach buildSparseReachND(const SystemND& sys, double prune) {
 
     SparseReach out; out.nCells = (int)N; out.nnz = 0;
     out.model.assign((size_t)N + 2, {}); out.targets.insert(TARGET);
+    out.actions = actions;
 
     auto cellLoDim = [&](int i, int j) { return sys.xlb[i] + j * sys.eta[i]; };
     auto isTargetMi = [&](const std::vector<int>& mi) {
@@ -188,9 +189,11 @@ SparseReach buildSparseReachND(const SystemND& sys, double prune) {
                 if (i == dx) break;
             }
 
-            double sumLo = 0, sumHi = 0;
-            for (const auto& iv : row) { sumLo += iv.lo; sumHi += iv.hi; }
-            row.push_back({SINK, std::max(0.0, 1.0 - sumHi), std::min(1.0, 1.0 - sumLo)});
+            // mass leaving the grid -> SINK, bounded tightly via the whole-grid-box
+            // complement (product of per-dim in-grid probabilities).
+            double gl = 1.0, gh = 1.0;
+            for (int i = 0; i < dx; ++i) { Bound g = transitionInterval1D(muLo[i], muHi[i], sys.sigma[i], sys.xlb[i], sys.xub[i]); gl *= g.lo; gh *= g.hi; }
+            row.push_back({SINK, std::max(0.0, 1.0 - gh), std::min(1.0, 1.0 - gl)});
             out.nnz += (long long)row.size();
             out.model[lin].push_back(std::move(row));
         }

src/abstraction.h

@@ -63,6 +63,9 @@ namespace abstraction {
         std::set<int> targets;      // { N }
         int nCells;
         long long nnz;              // total stored successor intervals (sparsity metric)
+        // Input vector for each action index (same action set for every non-target
+        // cell); lets a synthesized policy be simulated on the continuous system.
+        std::vector<std::vector<double>> actions;
     };
 
     // Build the sparse IMDP. `prune` drops successor cells whose upper bound <= prune