130 lines
6.7 KiB
Python
130 lines
6.7 KiB
Python
"""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()
|