diff --git a/tests/validation/multigrid/spacetime_union_plan.py b/tests/validation/multigrid/spacetime_union_plan.py index 223f857..81e1465 100644 --- a/tests/validation/multigrid/spacetime_union_plan.py +++ b/tests/validation/multigrid/spacetime_union_plan.py @@ -48,6 +48,20 @@ def _first_entry(covered, origin, child, c): return None +def _crossing_time(covered, origin, position, c, axes): + """Fig.11 convex T uses the last plane; Fig.13 concave T uses the first. + + Both use the same KKT formulation on different boundary geometry. + DOI/primary PDF: https://www.osti.gov/servlets/purl/1226215, §3.5.1. + """ + times=[(origin[j]+int(c[j]>0)-position[j])/c[j] for j in axes] + if len(axes)==1: return times[0] + sides=[bool(covered[tuple(origin+np.eye(2,dtype=int)[j]*c[j])]) for j in axes] + if sides[0]!=sides[1]: + raise UnresolvedTopology('mixed_corner_support',source_yx=origin.tolist()) + return min(times) if all(sides) else max(times) + + class UnionPlan: """Two-level, level-wide square geometry; all coordinates are level-local. @@ -72,9 +86,10 @@ class UnionPlan: if not np.array_equal(active, pair.fine_kind == 0): raise UnresolvedTopology('non_pair_active_coverage') self.maps = SimpleNamespace(coarse_active=~covered, coarse_covered=covered, fine_active=active) - self.outgoing = [[], []] - # Check transit FIRST: a fine exit may not be silently routed into a - # covered coarse cell. Returning fine packets need a separately derived rule. + self.outgoing, transit = [[], []], [] + # Guzik et al. §3.5.4, Fig.19 (https://www.osti.gov/servlets/purl/1226215): + # retain substep-0 outgoing G at source, then fill the substep-1 ghost. + # Rows are (q, fine source, next ghost); no covered coarse receiver credit. for q in range(1, 9): c = np.array((ref.cpu.CY[q], ref.cpu.CX[q])) boundary = active & ~np.roll(active, tuple(-c), axis=(0, 1)) @@ -83,10 +98,14 @@ class UnionPlan: for m in (0, 1): target = (p+(2-m)*c)//2 if covered[tuple(target)]: - raise UnresolvedTopology('fine_exit_reentry', q=q, substep=m, - fine_path_yx=[(p+j*c).tolist() for j in range(3-m)], coarse_target_yx=target.tolist()) - self.outgoing[m].append((q, int(y*self.nf+x), int(target[0]*nx+target[1]))) + if m != 0 or not active[tuple(p+2*c)]: + raise UnresolvedTopology('unsupported_transit', q=q, substep=m, source_yx=p.tolist()) + ghost=p+c + transit.append((q,int(y*self.nf+x),int(ghost[0]*self.nf+ghost[1]))) + else: + self.outgoing[m].append((q, int(y*self.nf+x), int(target[0]*nx+target[1]))) self.outgoing = [np.array(x, dtype=np.int64).reshape(-1, 3) for x in self.outgoing] + self.transit = np.array(transit,dtype=np.int64).reshape(-1,3) self.groups = [] for iy, ix in zip(*np.nonzero(~covered)): origin = np.array((iy, ix)) @@ -112,9 +131,7 @@ class UnionPlan: points = [] for p, entry, event in zip(starts, entries, events): position = (p+.5)/2 - planes = origin+(c > 0) - times = [(planes[j]-position[j])/c[j] for j in axes] - time = max(times) + time = _crossing_time(covered,origin,position,c,axes) if entry is not None and (4*time != entry[0] or event is None): raise UnresolvedTopology('nonconvex_entry_support', source_yx=origin.tolist(), q=q, quarter_time=entry[0], support_quarter_time=float(4*time)) @@ -140,6 +157,14 @@ def validate_ownership(plan): """Independent mask-based event coverage and four-quarter receiver accounting.""" n, nf, active = plan.n, plan.nf, plan.maps.fine_active actual = [(m,g.q,slot) for g in plan.groups for e in g.events if e is not None for m,slot in [e]] + expected_transit={(q,int(y*nf+x),int((y+ref.cpu.CY[q])*nf+x+ref.cpu.CX[q])) + for q in range(1,9) for y,x in zip(*np.nonzero(active)) + if not active[y+ref.cpu.CY[q],x+ref.cpu.CX[q]] + and active[y+2*ref.cpu.CY[q],x+2*ref.cpu.CX[q]]} + routes=[tuple(map(int,row)) for row in plan.transit] + if len(routes)!=len(set(routes)) or set(routes)!=expected_transit: + raise ValueError('transit provenance/uniqueness/coverage failed') + actual += [(1,q,ghost) for q,source,ghost in routes] expected = {(m,q,int((y-ref.cpu.CY[q])*nf+x-ref.cpu.CX[q])) for m in (0,1) for q in range(1,9) for y,x in zip(*np.nonzero(active)) if not active[y-ref.cpu.CY[q],x-ref.cpu.CX[q]]} @@ -178,7 +203,7 @@ def validate_ownership(plan): if not np.array_equal(g.samples,indices): raise ValueError('sample anchor provenance failed') for k,(dy,dx) in enumerate(((0,0),(0,1),(1,0),(1,1))): position=origin+(np.array((dy,dx))+.5)/2 - time=max((origin[j]+int(c[j]>0)-position[j])/c[j] for j in axes) + time=_crossing_time(plan.maps.coarse_covered,origin,position,c,axes) offset=position+c*time-(origin+.5+c*.5) point=(float(offset[1]),float(offset[0]),float(time-.5)) if g.kind=='corner' else (float(offset[1-axes[0]]),float(time-.5)) if g.points[k]!=point: raise ValueError('evaluation point provenance failed') @@ -193,7 +218,8 @@ def validate_ownership(plan): for m,rows in enumerate(plan.outgoing): expected_out = {(q,int(y*nf+x),int(((y+(2-m)*ref.cpu.CY[q])//2)*n+(x+(2-m)*ref.cpu.CX[q])//2)) for q in range(1,9) for y,x in zip(*np.nonzero(active)) - if not active[y+ref.cpu.CY[q],x+ref.cpu.CX[q]]} + if not active[y+ref.cpu.CY[q],x+ref.cpu.CX[q]] + and plan.maps.coarse_active[(y+(2-m)*ref.cpu.CY[q])//2,(x+(2-m)*ref.cpu.CX[q])//2]} actual_out = [tuple(map(int,row)) for row in rows] if len(actual_out)!=len(set(actual_out)) or set(actual_out)!=expected_out: raise ValueError('outgoing uniqueness/coverage failed') @@ -207,6 +233,7 @@ def validate_ownership(plan): if any(len(row)!=4 for row in packets.values()) or len(tokens)!=len(set(tokens)): raise ValueError('receiver needs four distinct quarter-packets') return dict(topology='supported',incoming_events=len(actual),outgoing_events=sum(map(len,plan.outgoing)), + transit_events=len(routes), groups=len(plan.groups),receivers=len(packets),receiver_quarter_packets=len(tokens), interpolation_accuracy=plan.interpolation_accuracy,runtime_supported=False, scope='pair-local ownership/coefficients; complete parent domain assumed, nested coverage unvalidated') diff --git a/tests/validation/multigrid/test_spacetime_union_plan.py b/tests/validation/multigrid/test_spacetime_union_plan.py index 2aab543..3e16435 100644 --- a/tests/validation/multigrid/test_spacetime_union_plan.py +++ b/tests/validation/multigrid/test_spacetime_union_plan.py @@ -14,6 +14,7 @@ import spacetime_gpu as frozen # pack_plan only: no CUDA import/owner construct from spacetime_table_guard import exact_weights MATRIX=[] +ENUMERATION={} CASES={ 'rectangle': ([(4,4,12,12)],True), 'two_separated': ([(3,3,6,6),(9,9,12,12)],True), @@ -21,11 +22,11 @@ CASES={ 'one_cell_diagonal_gap': ([(3,3,6,6),(7,7,10,10)],True), 'staggered_gap': ([(3,3,6,8),(7,5,11,10)],True), 'thin_rectangle': ([(4,4,5,11)],True), - 'concave_L': ([(4,4,5,7),(4,5,7,7)],False), - 'corner_contact': ([(3,3,6,6),(6,6,9,9)],False), - 'cross': ([(3,6,13,9),(6,3,9,13)],False), - 'narrow_notch': ([(3,3,6,11),(7,3,10,11),(3,9,10,11)],False), - 'hole': ([(3,3,13,5),(3,11,13,13),(3,5,5,11),(11,5,13,11)],False), + 'concave_L': ([(4,4,5,7),(4,5,7,7)],True), + 'corner_contact': ([(3,3,6,6),(6,6,9,9)],True), + 'cross': ([(3,6,13,9),(6,3,9,13)],True), + 'narrow_notch': ([(3,3,6,11),(7,3,10,11),(3,9,10,11)],True), + 'hole': ([(3,3,13,5),(3,11,13,13),(3,5,5,11),(11,5,13,11)],True), } @@ -44,6 +45,9 @@ def assert_same(test,a,b): for x,y in zip(a.groups,b.groups): test.assertEqual((x.q,x.source,x.receiver,x.kind,x.points,x.events), (y.q,y.source,y.receiver,y.kind,y.points,y.events)) + # Frozen pack_plan has no transit column and is NOT a union runtime format. + np.testing.assert_array_equal(getattr(a,'transit',np.empty((0,3),int)), + getattr(b,'transit',np.empty((0,3),int))) class UnionContract(unittest.TestCase): @@ -63,6 +67,9 @@ class UnionContract(unittest.TestCase): assert_same(self,baseline,UnionPlan(pair(order))) two=CASES['two_separated'][0] assert_same(self,UnionPlan(pair(two)),UnionPlan(pair(list(reversed(two))))) + lshape=CASES['concave_L'][0] + for decomposition in (list(reversed(lshape)),[(4,4,5,5),(4,5,7,7)],lshape+lshape): + assert_same(self,UnionPlan(pair(lshape)),UnionPlan(pair(decomposition))) def test_supported_unresolved_matrix(self): for name,(rectangles,supported) in CASES.items(): @@ -74,25 +81,113 @@ class UnionContract(unittest.TestCase): self.assertEqual(plan.proof['receiver_quarter_packets'],4*len(plan.receivers)) else: with self.assertRaises(UnresolvedTopology) as caught: UnionPlan(geometry) - self.assertEqual(caught.exception.code,'fine_exit_reentry') row.update(topology='unresolved',code=caught.exception.code,witness=caught.exception.detail) MATRIX.append(row) def test_L_witness_and_reflections(self): rectangles=CASES['concave_L'][0] - with self.assertRaises(UnresolvedTopology) as caught: UnionPlan(pair(rectangles)) - self.assertEqual(caught.exception.detail,dict(q=5,substep=0, - fine_path_yx=[[8,9],[9,10],[10,11]],coarse_target_yx=[5,5])) + plan=UnionPlan(pair(rectangles)) + self.assertIn((5,8*32+9,9*32+10),map(tuple,plan.transit)) directions=set() for reflect in (False,True): transformed=[(16-x1,y0,16-x0,y1) if reflect else (x0,y0,x1,y1) for x0,y0,x1,y1 in rectangles] for _ in range(4): - with self.assertRaises(UnresolvedTopology) as caught: UnionPlan(pair(transformed)) - detail=caught.exception.detail; directions.add(detail['q']) - p=np.array(detail['fine_path_yx']); geometry=pair(transformed) - self.assertEqual([bool(geometry.fine_kind[tuple(z)]==0) for z in p],[True,False,True]) + plan=UnionPlan(pair(transformed)); self.assertEqual(len(plan.transit),2) + for q,source,ghost in plan.transit: + directions.add(q); p=np.array(divmod(source,32)); c=np.array((ref.cpu.CY[q],ref.cpu.CX[q])) + self.assertEqual(ghost,(p+c)[0]*32+(p+c)[1]) + self.assertEqual([bool(plan.maps.fine_active[tuple(p+j*c)]) for j in range(3)],[True,False,True]) transformed=[(16-y1,x0,16-y0,x1) for x0,y0,x1,y1 in transformed] - self.assertGreater(len(directions),1) + self.assertEqual(directions,{5,6,7,8}) + + def test_Fig19_two_hop_transport_and_quarter_budget(self): + # Guzik §3.5.4 Fig.19: https://www.osti.gov/servlets/purl/1226215 . + # Isolated tagged transport only: no collision law or union evolution. + for name in ('concave_L','corner_contact','hole'): + p=UnionPlan(pair(CASES[name][0])); q,source,ghost=p.transit.T + stride=ref.cpu.CY[q]*p.nf+ref.cpu.CX[q]; target=ghost+stride + tags=(np.arange(len(q),dtype=np.float32)+1)/256 + active=p.maps.fine_active + def run(shifted,mutation=None): + scratch=np.full((9,p.nf,p.nf),np.nan,np.float32) + scratch.reshape(9,-1)[q,source]=tags + def moment(field): + value=field.astype(float)+(ref.cpu.W[:,None,None] if shifted else 0) + return .25*np.where(np.isfinite(field)&active,value,0).sum(axis=(1,2)) + def pull(field): + moved=np.stack([np.roll(field[j],(ref.cpu.CY[j],ref.cpu.CX[j]),axis=(0,1)) for j in range(9)]) + moved[:,~active]=np.nan + return moved + initial=moment(scratch) + carry=[(int(j),int(h),v.copy()) for j,h,v in zip(q,ghost,scratch.reshape(9,-1)[q,source])] + if mutation=='drop': carry.pop() + if mutation=='duplicate': carry.append(carry[0]) + first=pull(scratch) # m0 ghost arrivals are outside the ACTIVE mask. + inventory=np.zeros(9) + for j,h,value in carry: inventory[j]+=.25*(float(value)+(ref.cpu.W[j] if shifted else 0)) + np.testing.assert_allclose(moment(first)+inventory,initial,rtol=0,atol=2e-14) + if mutation=='early': + for j,h,value in carry: scratch.reshape(9,-1)[j,h]=value + scratch[:]=np.nan # m1 collision scratch reuse destroys an early injection. + while carry: + j,h,value=carry.pop() + if mutation!='early': scratch.reshape(9,-1)[j,h]=value + final=pull(scratch) + self.assertEqual(carry,[]) # Consumed route inventory, not subtraction of identical sums. + np.testing.assert_allclose(moment(final),initial,rtol=0,atol=2e-14) + arrivals=final.reshape(9,-1)[q,target] + np.testing.assert_array_equal(arrivals.view(np.uint32),tags.view(np.uint32)) + for shifted in (False,True): + run(shifted) + for mutation in ('drop','duplicate','early'): + with self.assertRaises(AssertionError): run(shifted,mutation) + for direction,slot,_ in p.transit: + self.assertNotIn((direction,slot),{tuple(row[:2]) for row in p.outgoing[0]}) + for change in (lambda b:setattr(b,'transit',b.transit[1:]), + lambda b:setattr(b,'transit',np.vstack((b.transit,b.transit[0]))), + lambda b:b.transit.__setitem__((0,2),b.transit[0,2]+1), + lambda b:b.transit.__setitem__((0,0),1)): + bad=copy.deepcopy(p); change(bad) + with self.assertRaises(ValueError): validate_ownership(bad) + bad=copy.deepcopy(p); direction,slot,_=p.transit[0] + y,x=divmod(slot,p.nf); cy,cx=ref.cpu.CY[direction],ref.cpu.CX[direction] + covered_target=((y+2*cy)//2)*p.n+(x+2*cx)//2 + bad.outgoing[0]=np.vstack((bad.outgoing[0],[direction,slot,covered_target])) + with self.assertRaises(ValueError): validate_ownership(bad) + + def test_Fig13_first_entry_and_wrong_crossing(self): + p=UnionPlan(pair(CASES['concave_L'][0])) + g=next(g for g in p.groups if g.source==4*16+5 and g.q==8) + self.assertEqual(g.kind,'corner') + self.assertEqual(sorted(point[-1] for point in g.points),[-.25,-.25,-.25,.25]) + # For a concave T first entry is min, not the frozen convex max. + damaged=copy.deepcopy(p); wrong=next(x for x in damaged.groups if x.source==g.source and x.q==g.q) + origin=np.array((4,5)); c=np.array((ref.cpu.CY[g.q],ref.cpu.CX[g.q])); points=[] + for child in ((0,0),(0,1),(1,0),(1,1)): + position=origin+(np.array(child)+.5)/2 + time=max((origin+ (c>0)-position)/c); offset=position+c*time-(origin+.5+c*.5) + points.append((float(offset[1]),float(offset[0]),float(time-.5))) + wrong.points=tuple(points) + with self.assertRaises(ValueError): validate_ownership(damaged) + + def test_all_511_small_masks(self): + counts=dict(masks=0,with_transit=0,transit_paths=0) + for bits in range(1,512): + rectangles=[(4+x,4+y,5+x,5+y) for y in range(3) for x in range(3) if bits & (1<<(3*y+x))] + p=UnionPlan(pair(rectangles)); counts['masks']+=1 + counts['with_transit']+=bool(len(p.transit)); counts['transit_paths']+=len(p.transit) + # Independent first-entry check: fine-cell membership before/after + # each quarter crossing, without the plan's plane min/max helper. + for g in p.groups: + c=np.array((ref.cpu.CY[g.q],ref.cpu.CX[g.q])); origin=np.array(divmod(g.source,p.n)) + for k,child in enumerate(((0,0),(0,1),(1,0),(1,1))): + start=2*origin+child; hits=[] + for t4 in (1,3): + before=(4*start+2+c*(2*t4-1))//4 + after=(4*start+2+c*(2*t4+1))//4 + if not p.maps.fine_active[tuple(before)] and p.maps.fine_active[tuple(after)]: hits.append(t4/4) + if hits: self.assertEqual(g.points[k][-1]+.5,hits[0]) + self.assertEqual(counts,dict(masks=511,with_transit=384,transit_paths=2048)); ENUMERATION.update(counts) def test_ordered_provenance_and_packet_mutations_fail(self): baseline=UnionPlan(pair(CASES['two_separated'][0])) @@ -112,7 +207,7 @@ class UnionContract(unittest.TestCase): def test_coefficients_and_directional_transfer_budget(self): rng=np.random.default_rng(31) - for name in ('two_separated','one_cell_face_gap','thin_rectangle'): + for name in CASES: p=UnionPlan(pair(CASES[name][0])); total=np.zeros(9) receiver={tuple(row):[] for row in p.receivers}; first_leftover=None for g in p.groups: @@ -187,5 +282,6 @@ class UnionContract(unittest.TestCase): if __name__=='__main__': result=unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.loadTestsFromTestCase(UnionContract)) print(json.dumps(dict(passed=result.wasSuccessful(),tests=result.testsRun, - scope='host-only topology/provenance/coefficient checks; no runtime or accuracy admission',matrix=MATRIX),indent=2)) + scope='host-only topology/provenance/coefficient checks; no runtime or accuracy admission',matrix=MATRIX, + enumeration=ENUMERATION),indent=2)) sys.exit(not result.wasSuccessful())