"""CPU-only analysis of the existing 521 frames; no solver or CUDA imports.""" import csv import hashlib import json import os from pathlib import Path import resource import sys import h5py import numpy as np ROOT = Path('/home/frank14f/CelerisLab') OUT = ROOT/'research/j3-offline-2026-10-12' SOURCE = ROOT/'research/j3-lag-2026-10-12/trajectory.h5' EXPECTED_SHA = '2b51f8ec5e0cc8ca72c73ac1b51dd611158368e6d6de791d85150fa111ea1986' SCALE = .002 def digest(path): with Path(path).open('rb') as stream: return hashlib.file_digest(stream, 'sha256').hexdigest() if hasattr(hashlib, 'file_digest') else hashlib.sha256(stream.read()).hexdigest() def main(): assert Path.cwd() == ROOT assert sys.executable == '/home/frank14f/anaconda3/envs/pycuda_3_10/bin/python' assert digest(SOURCE) == EXPECTED_SHA assert not (OUT/'results.json').exists(), 'Do not overwrite prior results' with h5py.File(SOURCE, 'r') as hf: step = hf['step'][:] u, v = hf['u'][:], hf['v'][:] assert u.shape == v.shape == (521, 33, 48) assert np.array_equal(step, np.arange(20760, 21281)) even = step % 2 == 0 ne, no = int(even.sum()), int((~even).sum()) assert (ne, no) == (261, 260) y = np.arange(1, 34) regions = dict(all=y > 0, bottom=y == 1, top=y == 33, interior=(y >= 2) & (y <= 32), top4=y >= 30) report = dict(input=str(SOURCE), input_sha256=EXPECTED_SHA, script_sha256=digest(__file__), python=sys.executable, velocity_scale=SCALE, even_frames=ne, odd_frames=no, lags={}, regions=[], parity={}, solver_steps=0, j3_acceptance='failed_unchanged') maps = {} for name, field in (('u', u), ('v', v)): for lag in (2, 4, 8): difference = field[lag:]-field[:-lag] rms = np.sqrt(np.mean(difference*difference, axis=(1, 2))) entry = dict(origins=len(difference), rms_mean=float(rms.mean()), rms_min=float(rms.min()), rms_max=float(rms.max()), max_abs=float(np.abs(difference).max()), signed_mean=float(difference.mean())) entry['rms_mean_scaled'] = entry['rms_mean']/SCALE entry['max_abs_scaled'] = entry['max_abs']/SCALE report['lags'].setdefault(str(lag), {})[name] = entry mean = field.mean(axis=0) variance = np.mean((field-mean)**2, axis=0) delta = field[even].mean(axis=0)-field[~even].mean(axis=0) between = (ne/521)*(no/521)*delta**2 # Independent weighted form checks the unequal group sizes. weighted = (ne/521)*(field[even].mean(axis=0)-mean)**2 + (no/521)*(field[~even].mean(axis=0)-mean)**2 np.testing.assert_allclose(between, weighted, rtol=1e-10, atol=1e-28) maps[name] = (np.sqrt(variance), delta) for region, mask in regions.items(): energy = float(variance[mask].sum()) fraction = float(mask.mean()) d = delta[mask] delta_energy = float(np.mean(d*d)) row = dict(component=name, region=region, nodes=int(mask.sum()*48), node_fraction=fraction, rms=float(np.sqrt(variance[mask].mean())), energy_fraction=energy/float(variance.sum()), parity_explained=float(between[mask].sum()/energy), parity_delta_rms=float(np.sqrt(delta_energy)), parity_delta_max_abs=float(np.abs(d).max()), parity_delta_signed_mean=float(d.mean()), parity_uniform_fraction=float(d.mean()**2/delta_energy)) row['concentration'] = row['energy_fraction']/fraction report['regions'].append(row) report['parity'][name] = next(row for row in report['regions'] if row['component'] == name and row['region'] == 'all') print(json.dumps({k:report[k] for k in ('lags', 'regions')}), flush=True) (OUT/'results.json').write_text(json.dumps(report, indent=2)+'\n') with (OUT/'regions.csv').open('w', newline='') as stream: writer = csv.DictWriter(stream, fieldnames=list(report['regions'][0])) writer.writeheader() writer.writerows(report['regions']) lag_rows = [dict(lag=int(lag), component=name, **entry) for lag, components in report['lags'].items() for name, entry in components.items()] with (OUT/'lags.csv').open('w', newline='') as stream: writer = csv.DictWriter(stream, fieldnames=list(lag_rows[0])) writer.writeheader() writer.writerows(lag_rows) # Cache is confined to this output directory, never the global environment. assert os.environ.get('MPLCONFIGDIR') == str(OUT/'.plot-cache') import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 2, figsize=(9, 7), constrained_layout=True) rms_limit = max(maps[name][0].max() for name in maps)/1e-7 parity_limit = max(np.abs(maps[name][1]).max() for name in maps)/1e-7 for col, name in enumerate(('u', 'v')): for row in (0, 1): ax = axes[row, col] image = ax.imshow(maps[name][row]/1e-7, origin='lower', extent=(-.5, 47.5, .5, 33.5), interpolation='nearest', aspect='equal', cmap='viridis' if row == 0 else 'RdBu_r', vmin=0 if row == 0 else -parity_limit, vmax=rms_limit if row == 0 else parity_limit) ax.set_title(f'{name}: temporal RMS' if row == 0 else f'{name}: even mean - odd mean') ax.axhline(29.5, color='black', linestyle='--', linewidth=.7) ax.plot([15.5, 31.5], [1, 1], color='black', linewidth=2) ax.set(xlabel='x (lattice units)', ylabel='y (lattice units)', yticks=[1, 16, 30, 33], yticklabels=['1 bottom', '16', '30', '33 top']) fig.colorbar(axes[0, 0].images[0], ax=axes[0, :], label='velocity / 1e-7', shrink=.85) fig.colorbar(axes[1, 0].images[0], ax=axes[1, :], label='signed velocity difference / 1e-7', shrink=.85) fig.suptitle('Saved J3 trajectory: steps 20760-21280 (521 frames)\nDashed: top four rows; black bottom segment: slot', fontsize=12) fig.savefig(OUT/'localization.png', dpi=160) plt.close(fig) report['peak_rss_mib'] = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024 report['figure_sha256'] = digest(OUT/'localization.png') assert digest(SOURCE) == EXPECTED_SHA assert not any(n == 'CelerisLab' or n.startswith('pycuda') or n == 'run_j1' for n in sys.modules) assert report['peak_rss_mib'] < 1024 (OUT/'results.json').write_text(json.dumps(report, indent=2)+'\n') print('CPU-only analysis complete; peak RSS MiB:', report['peak_rss_mib'], flush=True) if __name__ == '__main__': main()