@@ -74,19 +74,17 @@ SparseReach buildSparseReach1D(const System1D& sys, double prune) {
7474 double muLo = std::min (m1, m2), muHi = std::max (m1, m2);
7575
7676 solve::ActionDist row;
77- // mass into the target region [tlo,thi]
78- Bound bT = transitionInterval1D (muLo, muHi, sys.sigma , sys.tlo , sys.thi );
79- if (bT.hi > prune) row.push_back ({TARGET , bT.lo , bT.hi });
80-
81- // mass into non-target grid cells within the support window
77+ // Transitions to the (disjoint) grid cells in the support window; a target
78+ // cell routes its mass to the absorbing TARGET. No separate target-region
79+ // aggregate (which would double-count cells partially overlapping a
80+ // grid-unaligned target). See ISSUE-0008.
8281 int jmin = (int )std::floor ((muLo - W - sys.xlb ) / sys.eta );
8382 int jmax = (int )std::floor ((muHi + W - sys.xlb ) / sys.eta );
8483 jmin = std::max (0 , jmin); jmax = std::min (N - 1 , jmax);
8584 for (int j = jmin; j <= jmax; ++j) {
86- if (isTargetCell (j)) continue ; // counted in the TARGET aggregate
8785 double cl = cellLo (j), cr = cl + sys.eta ;
8886 Bound b = transitionInterval1D (muLo, muHi, sys.sigma , cl, cr);
89- if (b.hi > prune) row.push_back ({j, b.lo , b.hi });
87+ if (b.hi > prune) row.push_back ({isTargetCell (j) ? TARGET : j, b.lo , b.hi });
9088 }
9189
9290 // mass leaving the grid -> SINK (value 0), bounded TIGHTLY as the
@@ -104,22 +102,37 @@ SparseReach buildSparseReach1D(const System1D& sys, double prune) {
104102 return out;
105103}
106104
107- SparseReach buildSparseReachND (const SystemND& sys, double prune) {
108- const int dx = sys.dim_x , du = sys.dim_u ;
105+ // --- interval arithmetic (sound natural inclusion) for nonlinear mean bounds ----
106+ Ival operator +(const Ival& a, const Ival& b) { return Ival (a.lo + b.lo , a.hi + b.hi ); }
107+ Ival operator -(const Ival& a, const Ival& b) { return Ival (a.lo - b.hi , a.hi - b.lo ); }
108+ Ival operator *(const Ival& a, const Ival& b) {
109+ double p1 = a.lo *b.lo , p2 = a.lo *b.hi , p3 = a.hi *b.lo , p4 = a.hi *b.hi ;
110+ return Ival (std::min ({p1,p2,p3,p4}), std::max ({p1,p2,p3,p4}));
111+ }
112+ Ival operator +(const Ival& a, double s) { return Ival (a.lo + s, a.hi + s); }
113+ Ival operator *(double s, const Ival& a) { return s >= 0 ? Ival (s*a.lo , s*a.hi ) : Ival (s*a.hi , s*a.lo ); }
114+ Ival isquare (const Ival& a) {
115+ if (a.lo >= 0 ) return Ival (a.lo *a.lo , a.hi *a.hi );
116+ if (a.hi <= 0 ) return Ival (a.hi *a.hi , a.lo *a.lo );
117+ double m = std::max (-a.lo , a.hi );
118+ return Ival (0.0 , m*m);
119+ }
120+
121+ SparseReach buildSparseReachGeneral (const GridSpec& g, const MeanBoundFn& mean, double prune) {
122+ const int dx = g.dim_x , du = g.dim_u ;
109123 std::vector<int > Nd (dx);
110124 std::vector<long long > stride (dx);
111125 long long N = 1 ;
112126 for (int i = 0 ; i < dx; ++i) {
113- Nd[i] = std::max (1 , (int )std::llround ((sys .xub [i] - sys .xlb [i]) / sys .eta [i]));
127+ Nd[i] = std::max (1 , (int )std::llround ((g .xub [i] - g .xlb [i]) / g .eta [i]));
114128 stride[i] = N; N *= Nd[i];
115129 }
116130 const int TARGET = (int )N, SINK = (int )N + 1 ;
117131
118- // input actions = Cartesian product of per-dimension input grids
119132 std::vector<std::vector<double >> upts (du);
120133 for (int k = 0 ; k < du; ++k) {
121- int Mk = std::max (0 , (int )std::llround ((sys .uub [k] - sys .ulb [k]) / sys .ueta [k]));
122- for (int t = 0 ; t <= Mk; ++t) upts[k].push_back (sys .ulb [k] + t * sys .ueta [k]);
134+ int Mk = std::max (0 , (int )std::llround ((g .uub [k] - g .ulb [k]) / g .ueta [k]));
135+ for (int t = 0 ; t <= Mk; ++t) upts[k].push_back (g .ulb [k] + t * g .ueta [k]);
123136 }
124137 std::vector<std::vector<double >> actions;
125138 if (du == 0 ) actions.push_back ({});
@@ -138,61 +151,54 @@ SparseReach buildSparseReachND(const SystemND& sys, double prune) {
138151 out.model .assign ((size_t )N + 2 , {}); out.targets .insert (TARGET );
139152 out.actions = actions;
140153
141- auto cellLoDim = [&](int i, int j) { return sys .xlb [i] + j * sys .eta [i]; };
154+ auto cellLoDim = [&](int i, int j) { return g .xlb [i] + j * g .eta [i]; };
142155 auto isTargetMi = [&](const std::vector<int >& mi) {
143156 for (int i = 0 ; i < dx; ++i) {
144- double lo = cellLoDim (i, mi[i]), hi = lo + sys .eta [i];
145- if (!(lo >= sys .tlo [i] - 1e-12 && hi <= sys .thi [i] + 1e-12 )) return false ;
157+ double lo = cellLoDim (i, mi[i]), hi = lo + g .eta [i];
158+ if (!(lo >= g .tlo [i] - 1e-12 && hi <= g .thi [i] + 1e-12 )) return false ;
146159 }
147160 return true ;
148161 };
149162
150163 std::vector<int > mi (dx), wlo (dx), whi (dx), jt (dx);
151- std::vector<double > muLo (dx), muHi (dx);
164+ std::vector<double > muLo (dx), muHi (dx), cellLo (dx), cellHi (dx) ;
152165 for (long long lin = 0 ; lin < N; ++lin) {
153166 for (int i = 0 ; i < dx; ++i) mi[i] = (int )((lin / stride[i]) % Nd[i]);
154167 if (isTargetMi (mi)) { out.model [lin].push_back ({ {TARGET , 1.0 , 1.0 } }); out.nnz += 1 ; continue ; }
168+ for (int i = 0 ; i < dx; ++i) { cellLo[i] = cellLoDim (i, mi[i]); cellHi[i] = cellLo[i] + g.eta [i]; }
155169
156170 for (const auto & u : actions) {
157- for (int i = 0 ; i < dx; ++i) { // affine interval arithmetic for mean range
158- double lo = sys.c .empty () ? 0.0 : sys.c [i], hi = lo;
159- for (int j = 0 ; j < dx; ++j) {
160- double a = sys.A [i][j], xl = cellLoDim (j, mi[j]), xr = xl + sys.eta [j];
161- if (a >= 0 ) { lo += a * xl; hi += a * xr; } else { lo += a * xr; hi += a * xl; }
162- }
163- for (int k = 0 ; k < du; ++k) { lo += sys.B [i][k] * u[k]; hi += sys.B [i][k] * u[k]; }
164- muLo[i] = lo; muHi[i] = hi;
165- }
171+ mean (cellLo, cellHi, u, muLo, muHi); // SOUND per-dim mean enclosure
166172
173+ // Transitions to each grid cell in the per-dimension kernel window. Cells
174+ // are DISJOINT boxes, so there is no double counting. A window cell that is
175+ // a target cell routes its mass to the absorbing TARGET (no separate target
176+ // aggregate, which would double-count cells that partially overlap a
177+ // grid-unaligned target region).
167178 solve::ActionDist row;
168- { double tl = 1.0 , th = 1.0 ; // target box aggregate
169- for (int i = 0 ; i < dx; ++i) { Bound b = transitionInterval1D (muLo[i], muHi[i], sys.sigma [i], sys.tlo [i], sys.thi [i]); tl *= b.lo ; th *= b.hi ; }
170- if (th > prune) row.push_back ({TARGET , tl, th}); }
171-
172- bool any = true ; // per-dim kernel window
179+ bool any = true ;
173180 for (int i = 0 ; i < dx; ++i) {
174- double W = 6.0 * sys .sigma [i];
175- wlo[i] = std::max (0 , (int )std::floor ((muLo[i] - W - sys .xlb [i]) / sys .eta [i]));
176- whi[i] = std::min (Nd[i] - 1 , (int )std::floor ((muHi[i] + W - sys .xlb [i]) / sys .eta [i]));
181+ double W = 6.0 * g .sigma [i];
182+ wlo[i] = std::max (0 , (int )std::floor ((muLo[i] - W - g .xlb [i]) / g .eta [i]));
183+ whi[i] = std::min (Nd[i] - 1 , (int )std::floor ((muHi[i] + W - g .xlb [i]) / g .eta [i]));
177184 if (wlo[i] > whi[i]) any = false ;
178185 jt[i] = wlo[i];
179186 }
180- if (any) while (true ) { // Cartesian product over windows
181- if (!isTargetMi (jt)) {
182- double pl = 1.0 , ph = 1.0 ;
183- for (int i = 0 ; i < dx; ++i) { double cl = cellLoDim (i, jt[i]), cr = cl + sys.eta [i];
184- Bound b = transitionInterval1D (muLo[i], muHi[i], sys.sigma [i], cl, cr); pl *= b.lo ; ph *= b.hi ; }
185- if (ph > prune) { long long lj = 0 ; for (int i = 0 ; i < dx; ++i) lj += (long long )jt[i] * stride[i];
187+ if (any) while (true ) {
188+ double pl = 1.0 , ph = 1.0 ;
189+ for (int i = 0 ; i < dx; ++i) { double cl = cellLoDim (i, jt[i]), cr = cl + g.eta [i];
190+ Bound b = transitionInterval1D (muLo[i], muHi[i], g.sigma [i], cl, cr); pl *= b.lo ; ph *= b.hi ; }
191+ if (ph > prune) {
192+ if (isTargetMi (jt)) { row.push_back ({TARGET , pl, ph}); }
193+ else { long long lj = 0 ; for (int i = 0 ; i < dx; ++i) lj += (long long )jt[i] * stride[i];
186194 row.push_back ({(int )lj, pl, ph}); }
187195 }
188196 int i = 0 ; for (; i < dx; ++i) { if (++jt[i] <= whi[i]) break ; jt[i] = wlo[i]; }
189197 if (i == dx) break ;
190198 }
191199
192- // mass leaving the grid -> SINK, bounded tightly via the whole-grid-box
193- // complement (product of per-dim in-grid probabilities).
194- double gl = 1.0 , gh = 1.0 ;
195- 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 ; }
200+ double gl = 1.0 , gh = 1.0 ; // outside-grid via grid-box complement
201+ for (int i = 0 ; i < dx; ++i) { Bound gg = transitionInterval1D (muLo[i], muHi[i], g.sigma [i], g.xlb [i], g.xub [i]); gl *= gg.lo ; gh *= gg.hi ; }
196202 row.push_back ({SINK , std::max (0.0 , 1.0 - gh), std::min (1.0 , 1.0 - gl)});
197203 out.nnz += (long long )row.size ();
198204 out.model [lin].push_back (std::move (row));
@@ -203,5 +209,28 @@ SparseReach buildSparseReachND(const SystemND& sys, double prune) {
203209 return out;
204210}
205211
212+ SparseReach buildSparseReachND (const SystemND& sys, double prune) {
213+ GridSpec g;
214+ g.dim_x = sys.dim_x ; g.dim_u = sys.dim_u ;
215+ g.xlb = sys.xlb ; g.xub = sys.xub ; g.eta = sys.eta ;
216+ g.ulb = sys.ulb ; g.uub = sys.uub ; g.ueta = sys.ueta ;
217+ g.sigma = sys.sigma ; g.tlo = sys.tlo ; g.thi = sys.thi ;
218+ const auto A = sys.A ; const auto B = sys.B ; const auto c = sys.c ;
219+ const int dx = sys.dim_x , du = sys.dim_u ;
220+ MeanBoundFn affine = [A, B, c, dx, du](const std::vector<double >& cl, const std::vector<double >& ch,
221+ const std::vector<double >& u,
222+ std::vector<double >& muLo, std::vector<double >& muHi) {
223+ muLo.assign (dx, 0.0 ); muHi.assign (dx, 0.0 );
224+ for (int i = 0 ; i < dx; ++i) {
225+ double lo = c.empty () ? 0.0 : c[i], hi = lo;
226+ for (int j = 0 ; j < dx; ++j) { double a = A[i][j];
227+ if (a >= 0 ) { lo += a * cl[j]; hi += a * ch[j]; } else { lo += a * ch[j]; hi += a * cl[j]; } }
228+ for (int k = 0 ; k < du; ++k) { lo += B[i][k] * u[k]; hi += B[i][k] * u[k]; }
229+ muLo[i] = lo; muHi[i] = hi;
230+ }
231+ };
232+ return buildSparseReachGeneral (g, affine, prune);
233+ }
234+
206235} // namespace abstraction
207236} // namespace impact
0 commit comments