test: accumulate same-trajectory space-time GPU rounding ledger

tests/validation/multigrid/run_spacetime_gpu_n64.py

@@ -23,6 +23,7 @@ ref = gpu.ref
 CPU_SHA = '180eef90a1c3f43113fe18cc54981ee4723cfc0cc2724334907ff188022fa770'
 UNIFORM_SHA = '01e7415ae39456f2de7f0c3d331afb971c49d19353b43038650e32e80737f1e7'
 CUDA_SHA = 'c47869688f7ac2ec93411c0af233887fe8b77405838299e4ae52c8963fd97fe4'
+AUDIT_SHA = 'fba340b1533bec308b666ee84666df0ac9d08e2872be31ee379e2d1447565158'
 STEPS = (256, 512, 768, 1024)
 TIMES = (.12, .24, .36, .48)
 FACTORS = dict(velocity_l2=.8, pressure_l2=1., energy_relative=1.)
@@ -56,14 +57,82 @@ def total(state, maps):
     return ref.cpu.MOMENTS.T @ ref.directional(ref.cpu.State(*arrays), maps)
 
 
-def warm_and_reset(run, initial, check_cap):
+def differing_fields(actual, expected):
+    """Compare serialized sample numerics exactly, including the sign of zero."""
+    if actual.keys() != expected.keys():
+        return ['field_set']
+    return [key for key in actual if np.asarray(actual[key]).dtype != np.asarray(expected[key]).dtype
+            or np.asarray(actual[key]).shape != np.asarray(expected[key]).shape
+            or np.asarray(actual[key]).tobytes() != np.asarray(expected[key]).tobytes()]
+
+
+class LedgerTotals:
+    """Compensated directional sums; retain absolute step sums to expose cancellation."""
+    names = ('net_change','physical_collision','coarse_collision','fine0_collision','fine1_collision',
+             'child_rounding','child_interpolation','child_cast','register_rounding','receiver_rounding','unexplained')
+
+    def __init__(self):
+        self.sums = np.zeros((len(self.names),9)); self.correction = self.sums.copy()
+        self.absolute = self.sums.copy(); self.absolute_moments = np.zeros((len(self.names),3))
+        self.maxima = self.sums.copy(); self.max_steps = np.zeros(self.sums.shape,np.int64)
+        self.cells = np.zeros((3,4)); self.cell_correction = self.cells.copy()
+        self.steps = 0; self.bound = 0.; self.maximum = 0.
+
+    def add(self, ledger):
+        values = np.array([ledger[name] for name in self.names],np.float64)
+        cells = np.asarray(ledger['collision_cell_mass'],np.float64)
+        if (values.shape != self.sums.shape or not np.isfinite(values).all()
+                or cells.shape != (3,4) or not np.isfinite(cells).all()
+                or not np.isfinite(ledger['unexplained_bound']) or ledger['unexplained_bound'] <= 0):
+            raise AssertionError('invalid diagnostic ledger')
+        delta = values-self.correction
+        updated = self.sums+delta
+        self.correction = (updated-self.sums)-delta
+        self.sums = updated
+        self.absolute += abs(values)
+        self.absolute_moments += abs(values @ ref.cpu.MOMENTS)
+        larger = abs(values) > self.maxima
+        self.max_steps[larger] = self.steps+1; self.maxima = np.maximum(self.maxima,abs(values))
+        delta = cells-self.cell_correction
+        updated = self.cells+delta
+        self.cell_correction = (updated-self.cells)-delta; self.cells = updated
+        self.steps += 1; self.bound += ledger['unexplained_bound']
+        self.maximum = max(self.maximum,float(abs(values[-1]).max()))
+
+    def snapshot(self, endpoint_change=None):
+        values = dict(zip(self.names,self.sums))
+        accounted = sum(values[key] for key in ('physical_collision','child_rounding','register_rounding','receiver_rounding'))
+        closure = None if endpoint_change is None else endpoint_change-accounted
+        telescoping = None if endpoint_change is None else endpoint_change-values['net_change']
+        closure_pass = None if closure is None else max(float(abs(closure).max()),float(abs(telescoping).max())) <= self.bound
+        collision_split = values['physical_collision']-sum(values[key] for key in ('coarse_collision','fine0_collision','fine1_collision'))
+        child_split = values['child_rounding']-values['child_interpolation']-values['child_cast']
+        split_pass = max(float(abs(collision_split).max()),float(abs(child_split).max())) <= self.bound
+        return dict(intervals=self.steps, summed_per_step_bound=self.bound, max_step_unexplained=self.maximum,
+            closure_pass=closure_pass,component_split_pass=split_pass,
+            endpoint_status='not_measured' if endpoint_change is None else 'measured',
+            endpoint_directional_change=None if endpoint_change is None else endpoint_change.tolist(),
+            closure_residual=None if closure is None else closure.tolist(),
+            telescoping_residual=None if telescoping is None else telescoping.tolist(),
+            closure_mass_momentum=None if closure is None else (ref.cpu.MOMENTS.T@closure).tolist(),
+            collision_split_residual=collision_split.tolist(),child_split_residual=child_split.tolist(),
+            collision_cell_mass=dict(levels=['coarse','fine0','fine1'],columns=['net','positive','negative','l1'],
+                                     applied_volume_weights=[1.,.25,.25],values=self.cells.tolist()),
+            components={key:dict(directional=self.sums[i].tolist(),mass_momentum=(self.sums[i]@ref.cpu.MOMENTS).tolist(),
+                absolute_step_directional_sum=self.absolute[i].tolist(),
+                absolute_step_mass_momentum_sum=self.absolute_moments[i].tolist(),
+                max_absolute_interval=self.maxima[i].tolist(),max_absolute_interval_step=self.max_steps[i].tolist())
+                for i,key in enumerate(self.names)})
+
+
+def warm_and_reset(run, initial, check_cap, *, audit=False):
     """Time 32 synchronized intervals, then restore exact uploaded arrival F."""
     run.initialize(initial)
     before = run.download()
     started = time.monotonic()
     for _ in range(32):
         check_cap()
-        run.step(audit=False)
+        run.step(audit=audit)
     seconds = time.monotonic()-started
     warm_health = run.health_counts.tolist()
     run.initialize(initial)
@@ -73,7 +142,7 @@ def warm_and_reset(run, initial, check_cap):
     if run.step_count != 0 or after.tick != 0 or run.health_counts.any():
         raise AssertionError('warm-up reset did not restore tick/counters')
     return after, dict(intervals=32, seconds=seconds, estimated_1024_seconds=32*seconds,
-                       health_counts=warm_health, reset_verified=True)
+                       health_counts=warm_health, reset_verified=True, ledger_audit=audit)
 
 
 def assess(rows, cpu, uniform, health_counts):
@@ -114,6 +183,7 @@ def main():
     parser.add_argument('--gpu-admitted', action='store_true')
     parser.add_argument('--cpu-receipt', type=Path, required=True)
     parser.add_argument('--uniform-receipt', type=Path, required=True)
+    parser.add_argument('--audit-against', type=Path, help='audit every interval and compare exact saved 1779cf2 samples')
     parser.add_argument('--output', type=Path, required=True)
     args = parser.parse_args()
     if not args.gpu_admitted or not os.environ.get('CUDA_VISIBLE_DEVICES'):
@@ -121,6 +191,13 @@ def main():
     if args.output.exists():
         parser.error('output exists; preserve evidence and choose a new path')
     cpu, uniform = load_references(args.cpu_receipt, args.uniform_receipt)
+    previous = None
+    if args.audit_against is not None:
+        if sha(args.audit_against) != AUDIT_SHA:
+            raise ValueError('pinned GPU trajectory receipt differs')
+        previous = json.loads(args.audit_against.read_text())
+        if previous['source']['numpy'] != np.__version__:
+            raise RuntimeError('audited repeat requires unchanged NumPy environment')
     gpu.check_sources()
     if sha(Path(gpu.__file__).with_suffix('.cu')) != CUDA_SHA:
         raise RuntimeError('reviewed CUDA source changed')
@@ -145,10 +222,15 @@ def main():
             uniform_storage='historical_shifted_FP32'),
         frozen_screen=dict(campaign_limits=LIMITS, experimental_factors=FACTORS,
                            gpu_cpu_equivalence_tolerance=None, historical_pressure_threshold=.024506347612699403),
-        boundary_budget_status='not_measured', max_directional_residual=None,
+        boundary_budget_status='pending_measurement' if previous else 'not_measured', max_directional_residual=None,
         minima_scope='initial and four sampled physical arrival states; boundary minima at four samples',
         limits=dict(seconds=180., peak_rss_mib=512.), cuda_visible_devices=os.environ['CUDA_VISIBLE_DEVICES'])
     run = context = None
+    ledger_totals = LedgerTotals() if previous else None
+    if previous:
+        report['diagnostic'] = dict(completed=False, passed=False, reference_sha256=AUDIT_SHA,
+            reference_source=previous['source']['commit'], sample_checks=[], prefixes=[],
+            scope='net drift attribution; exact saved diagnostics only, not full-field trajectory certification')
     phase = 'setup'
 
     def save():
@@ -159,6 +241,11 @@ def main():
         if run is not None:
             report['health_counts'] = run.health_counts.tolist()
             report['completed_intervals'] = run.step_count
+        if ledger_totals is not None:
+            report['audited_intervals'] = ledger_totals.steps
+            if ledger_totals.steps:
+                report['boundary_budget_status'] = 'measured' if ledger_totals.steps==1024 else 'measured_prefix'
+                report['diagnostic']['accumulated'] = ledger_totals.snapshot()
         ref.write_receipt(args.output, report)
 
     def check_cap():
@@ -174,9 +261,13 @@ def main():
         report.update(device=context.name, sm_arch=context.sm_arch, build_metadata=run.build_metadata,
                       table_sha256=gpu.table_hash(run.tables))
         report['table_validation'] = table_guard.validate_tables(run.plan,run.tables)
+        if previous:
+            for key in ('device','sm_arch','table_sha256','build_metadata'):
+                if report[key] != previous[key]:
+                    raise AssertionError(f'audited repeat identity differs: {key}')
         initial = ref.cpu.State(tg.initial((64,64),0,64,.0192), tg.initial((128,128),1,64,.0192))
         phase = 'warm_reset'
-        initial, report['warmup'] = warm_and_reset(run, initial, check_cap)
+        initial, report['warmup'] = warm_and_reset(run, initial, check_cap, audit=previous is not None)
         maps = run.plan.maps
         masks = (maps.coarse_active, maps.fine_active)
         arrays = (initial.coarse_f, initial.fine_f)
@@ -184,13 +275,19 @@ def main():
         moment0 = total(initial, maps)
         report['initial'] = dict(mass=mass0, momentum=moment0[1:].tolist(),
             minimum_population=min(float(a[:,m].min()) for a,m in zip(arrays,masks)))
+        if previous and differing_fields(report['initial'],previous['initial']):
+            raise AssertionError('audited repeat initial state diagnostics differ')
+        direction0 = ref.directional(ref.cpu.State(*(a.astype(np.float64) for a in arrays)),maps)
         save()
         check_cap()
         phase = 'trajectory'
         trajectory_started = time.monotonic()
         for step in range(1,1025):
             check_cap()
-            run.step(audit=False)
+            ledger = run.step(audit=previous is not None)
+            if ledger_totals is not None:
+                ledger_totals.add(ledger)
+                report['max_directional_residual'] = ledger_totals.maximum
             if step in STEPS:
                 state = run.download()
                 if state.tick != 2*step:
@@ -207,6 +304,16 @@ def main():
                     report['failed_sample'] = {key: repr(value) for key,value in row.items()}
                     raise FloatingPointError('nonfinite sample diagnostic')
                 report['rows'].append(row)
+                if previous:
+                    changed = differing_fields(row,previous['rows'][len(report['rows'])-1])
+                    report['diagnostic']['sample_checks'].append(dict(step=step,differing_fields=changed))
+                    endpoint = ref.directional(ref.cpu.State(*(a.astype(np.float64) for a in arrays)),maps)-direction0
+                    prefix = ledger_totals.snapshot(endpoint)
+                    report['diagnostic']['prefixes'].append(prefix)
+                    if not prefix['closure_pass'] or not prefix['component_split_pass']:
+                        raise AssertionError('cumulative ledger does not close within summed existing per-step bounds')
+                    if changed:
+                        raise AssertionError(f'audited repeat sample differs: {changed}')
                 save()
                 print(json.dumps(row, allow_nan=False), flush=True)
         report['trajectory_seconds'] = time.monotonic()-trajectory_started
@@ -214,6 +321,8 @@ def main():
         report.update(assess(report['rows'], cpu, uniform, run.health_counts))
         report['completed'] = True
         report['passed'] = report['campaign_pass'] and report['experimental_pass']
+        if previous:
+            report['diagnostic'].update(completed=True,passed=ledger_totals.steps==1024 and len(report['rows'])==4)
     except Exception as exc:
         report['failure'] = dict(phase=phase, error=repr(exc))
         raise

tests/validation/multigrid/spacetime_gpu.py

@@ -83,13 +83,21 @@ def rounding_ledger(plan, tables, before, after, g0, predicted, fine_phases, pre
     """Audit actual stored FP32 phases; never call a missing packet 'rounding'."""
     maps, n = plan.maps, plan.n
     f64 = lambda a: np.asarray(a, dtype=np.float64)
-    collision = (f64(g0[:, maps.coarse_active])-f64(before.coarse_f[:, maps.coarse_active])).sum(axis=1)
+    difference = f64(g0[:, maps.coarse_active])-f64(before.coarse_f[:, maps.coarse_active])
+    coarse_collision = difference.sum(axis=1)
+    cell_mass = [difference.sum(axis=0)]
+    collision = coarse_collision.copy(); fine_collision = []
     for f, g in fine_phases:
-        collision += .25*(f64(g[:, maps.fine_active])-f64(f[:, maps.fine_active])).sum(axis=1)
+        difference = f64(g[:, maps.fine_active])-f64(f[:, maps.fine_active])
+        fine_collision.append(.25*difference.sum(axis=1))
+        cell_mass.append(.25*difference.sum(axis=0))
+        collision += fine_collision[-1]
     np.testing.assert_array_equal(h, pre.astype(np.float32))
     origins = f64(g0).ravel()[tables['origins']]
     np.testing.assert_allclose(pre.mean(axis=1), origins, rtol=0, atol=2e-14)
     group = np.bincount(tables['origins']//(n*n), weights=f64(h).mean(axis=1)-origins, minlength=9)
+    interpolation = np.bincount(tables['origins']//(n*n), weights=pre.mean(axis=1)-origins, minlength=9)
+    child_cast = np.bincount(tables['origins']//(n*n), weights=f64(h).mean(axis=1)-pre.mean(axis=1), minlength=9)
     packet = np.empty(tables['times'].shape, np.float64)
     for when, source in ((-1, h), (0, fine_phases[0][1]), (1, fine_phases[1][1])):
         mask = tables['times'] == when
@@ -107,11 +115,17 @@ def rounding_ledger(plan, tables, before, after, g0, predicted, fine_phases, pre
     register_error = np.bincount(q, weights=reg_error, minlength=9)
     # directional() must sum converted arrays, not accumulate in FP32.
     total = lambda s: ref.directional(ref.cpu.State(f64(s.coarse_f), f64(s.fine_f)), maps)
-    raw = total(after)-total(before)-collision
+    net = total(after)-total(before)
+    raw = net-collision
     unexplained = raw-group-store_error-register_error
     limit = float(64*np.finfo(float).eps*total(before).sum())
     np.testing.assert_allclose(unexplained, 0, rtol=0, atol=limit)
-    return dict(raw=raw.tolist(), physical_collision=collision.tolist(),
+    return dict(raw=raw.tolist(), net_change=net.tolist(), physical_collision=collision.tolist(),
+                coarse_collision=coarse_collision.tolist(), fine0_collision=fine_collision[0].tolist(),
+                fine1_collision=fine_collision[1].tolist(), child_interpolation=interpolation.tolist(),
+                child_cast=child_cast.tolist(),
+                collision_cell_mass=[[float(a.sum()),float(a[a>0].sum()),float(a[a<0].sum()),float(abs(a).sum())]
+                                     for a in cell_mass],
                 child_rounding=group.tolist(), receiver_rounding=store_error.tolist(),
                 register_rounding=register_error.tolist(), unexplained=unexplained.tolist(),
                 unexplained_bound=limit, minimum_pre=float(pre.min()), minimum_stored=float(h.min()))

tests/validation/multigrid/test_spacetime_gpu.py

@@ -102,6 +102,10 @@ class Host(unittest.TestCase):
         args = (p,t,state,state,g0,g0,phases,pre,h)
         ledger = gpu.rounding_ledger(*args,registers)
         self.assertEqual(max(map(abs,ledger['unexplained'])),0)
+        np.testing.assert_array_equal(ledger['physical_collision'],np.array(ledger['coarse_collision'])
+                                      +ledger['fine0_collision']+ledger['fine1_collision'])
+        np.testing.assert_array_equal(ledger['child_rounding'],np.array(ledger['child_cast'])+ledger['child_interpolation'])
+        np.testing.assert_array_equal(ledger['collision_cell_mass'],np.zeros((3,4)))
         registers[0] = .01
         with self.assertRaisesRegex(AssertionError,'register mismatch'):
             gpu.rounding_ledger(*args,registers)
@@ -149,6 +153,7 @@ class Host(unittest.TestCase):
         self.assertEqual(run.health_counts.tolist(),[0,0,0,0,1,0])
 
     def test_warmup_resets_arrivals_and_counters(self):
+        audits=[]
         class Fake:
             def initialize(self,state):
                 self.state=ref.cpu.State(state.coarse_f.copy(),state.fine_f.copy())
@@ -156,13 +161,40 @@ class Host(unittest.TestCase):
             def download(self):
                 return ref.cpu.State(self.state.coarse_f.copy(),self.state.fine_f.copy(),tick=2*self.step_count)
             def step(self,*,audit):
-                assert audit is False
+                audits.append(audit)
                 self.step_count+=1; self.state.coarse_f[:]+=1e-4
         fake=Fake(); check=Mock(); state=initial(self.plan)
         actual,warm=trajectory.warm_and_reset(fake,state,check)
         np.testing.assert_array_equal(actual.coarse_f,state.coarse_f)
         self.assertEqual(check.call_count,32); self.assertEqual(fake.step_count,0)
         self.assertTrue(warm['reset_verified']); self.assertEqual(warm['intervals'],32)
+        trajectory.warm_and_reset(fake,state,check,audit=True)
+        self.assertEqual(audits,[False]*32+[True]*32)
+
+    def test_accumulated_drift_ledger_and_exact_samples(self):
+        totals=trajectory.LedgerTotals()
+        for sign in (1.,-1.):
+            ledger={name:np.zeros(9) for name in totals.names}
+            ledger['coarse_collision'][0]=sign*2.**-20
+            ledger['physical_collision']=ledger['coarse_collision'].copy()
+            ledger['child_cast'][5]=sign*2.**-30
+            ledger['child_interpolation'][5]=sign*2.**-50
+            ledger['child_rounding']=ledger['child_cast']+ledger['child_interpolation']
+            ledger['receiver_rounding'][5]=-sign*2.**-32
+            ledger['net_change']=sum(ledger[k] for k in ('physical_collision','child_rounding','receiver_rounding'))
+            mass=ledger['coarse_collision'].sum()
+            ledger.update(unexplained_bound=1e-12,collision_cell_mass=[[mass,max(mass,0.),min(mass,0.),abs(mass)],[0.]*4,[0.]*4])
+            totals.add(ledger)
+        result=totals.snapshot(np.zeros(9))
+        self.assertEqual(result['intervals'],2)
+        self.assertEqual(result['components']['physical_collision']['directional'],[0.]*9)
+        self.assertGreater(result['components']['physical_collision']['absolute_step_mass_momentum_sum'][0],0)
+        self.assertEqual(result['components']['physical_collision']['max_absolute_interval_step'][0],1)
+        self.assertEqual(totals.snapshot()['endpoint_status'],'not_measured')
+        self.assertFalse(totals.snapshot(np.ones(9)*1e-8)['closure_pass'])
+        self.assertEqual(trajectory.differing_fields(dict(a=1.,b=[0.,2.]),dict(a=1.,b=[0.,2.])),[])
+        self.assertEqual(trajectory.differing_fields(dict(a=0.),dict(a=-0.)),['a'])
+        self.assertEqual(trajectory.differing_fields(dict(a=1.),dict(a=np.nextafter(1.,2.))),['a'])
 
     def test_trajectory_screens_differences_and_json(self):
         rows=[dict(tstar=t,velocity_l2=np.float64(.01),pressure_l2=.03,energy_relative=.01,