"""Offline Fig2 estimands; preserve raw same-clock errors beside registration. No thresholds are fitted here. Probe DTW is the actual Legacy helper, not a field metric; independent wake phase is a diagnostic, not a time correction. """ import argparse import gc import importlib.util import json from pathlib import Path import sys import types import numpy as np ROOT=Path(__file__).resolve().parents[2] sys.path.insert(0,str(ROOT)) from research.fig2_collection import run from research.fig2_materials.metrics import field_error, motion_summary def frequency(signal,dt): """Largest non-DC rectangular-window FFT bin, descriptive resolution only.""" a=np.asarray(signal,float);a=a-a.mean() if len(a)<3 or np.max(np.abs(a))<1e-12:return None power=np.abs(np.fft.rfft(a))**2 i=1+int(np.argmax(power[1:])) return dict(f_D_U=float(np.fft.rfftfreq(len(a),dt)[i]),bin_width_D_U=1/(len(a)*dt), method='demeaned rectangular window; dominant non-DC bin; not a confidence interval') def native_summary(rows,scene): start=max(0,len(rows)-640000) if scene!='taylor' else 0 a=rows[start:] y=a[:,1]/40.;cy=a[:,3:9].sum(axis=1)/(.5*.005**2*40) return dict(window_D_U=[start*.005/40,len(rows)*.005/40], motion=motion_summary(y),motion_frequency=frequency(y,.005/40), structural_RHS_Cy_RMS=float(np.sqrt(np.mean(cy**2))), structural_RHS_Cy_mean=float(cy.mean()), max_native_delta_Cy=float(np.abs(np.diff(cy)).max()) if len(cy)>1 else None, force_frequency=frequency(cy,.005/40), no_donor=float(rows[:,13].sum()),eq_only=float(rows[:,12].sum()), crossings=float(rows[:,14].sum()), force_convention='body force, sum pinball CME_y+Wen_y; native previous-step RHS, no filtering') def compare_selected(target,current,u0=.005): """Fixed lab ROI; all actual compared roles must be fluid and finite here.""" ys,xs=slice(311,711),slice(1360,2161) t,c=target['ux'][:,ys,xs],current['ux'][:,ys,xs] tv,cv=target['uy'][:,ys,xs],current['uy'][:,ys,xs] if t.shape!=c.shape or tv.shape!=cv.shape:raise ValueError('selected shape mismatch') for fields in (target,current): if not np.all((fields['flags'][:,ys,xs]&1)!=0):raise ValueError('nonfluid selected ROI') if not all(np.isfinite(a).all() for a in (t,c,tv,cv)):raise ValueError('nonfinite selected ROI') e2=((c.astype(float)-t)/u0)**2+((cv.astype(float)-tv)/u0)**2 return dict(per_slot_E=np.sqrt(e2.mean(axis=(1,2))),aggregate_E=float(np.sqrt(e2.mean()))) def main(): ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('directory',type=Path) args=ap.parse_args();p=args.directory case=json.loads((p/'case.json').read_text());term=json.loads((p/'terminal.json').read_text()) if case['plan']['smoke'] or term['status']!='three_roles_completed_not_qualified': raise ValueError('physical analysis requires all three complete roles; smoke/partial is unscored') helpers=run.load_helpers() # Reuse the canonical pure module under the same isolated read-only package. native=run.load_file('_fig2_readonly_drl.legacy_test.metrics', run.ACQUISITION.parent/'legacy_test/metrics.py') reports={k:json.loads((p/k/'report.json').read_text()) for k in ('target','zero','controlled')} sig={} result=dict(scene=case['scene'],source=case['source'],Re_D=50.,roles={}, primary='unaligned same-clock/event U-normalized vector errors; common lab ROI', periodic_phase='independently registered centre-uy phase diagnostic, not a shared-clock comparison', notes=['Modern moving-plant data; no Legacy numerical equivalence or physical qualification', 'Motion, native structural force, policy SI means and probe DTW are different estimands']) for role in ('target','zero','controlled'): with np.load(p/role/'signals.npz') as z:sig[role]={k:z[k] for k in z.files} if role!='target': with np.load(p/role/'native.npz') as z: result['roles'][role]=native_summary(z['rows'],case['scene']) t=sig['target']['sensors'] for role in ('zero','controlled'): s=sig[role]['sensors'] if not np.array_equal(sig[role]['absolute_steps'],sig['target']['absolute_steps']): raise ValueError('probe absolute clocks mismatch') result['roles'][role]['unaligned_probe_area_velocity_RMS_over_U']=float(np.sqrt(np.mean(((s/156-t/156)/.005)**2))) result['roles'][role]['probe_channel_correlations']=native.channel_correlations(t,s) if case['scene']=='taylor': spec=types.SimpleNamespace(scene='vortex',conv_len=30) dtw=native.native_similarity(spec,t,s,len(s)-1) result['roles'][role]['native_probe_DTW']=dict(similarity=dtw,definition='source vortex training-window roll; supplementary, not event field recovery',conv_len=30) own=helpers.select_vortex_event(s,'vortex_taylor') target_event=helpers.select_vortex_event(t,'vortex_taylor') result['roles'][role]['own_event_peak']=int(own['peak_index']) result['roles'][role]['event_peak_delay_D_U']=(int(own['peak_index'])-int(target_event['peak_index']))*.4 else: result['roles'][role]['native_probe_DTW']=native.cycle_similarity(t,s,30) result['roles'][role]['native_probe_DTW'].update(definition='canonical six-channel legacy cycle DTW, training D20 units', lag_column=1,target_template_indices=[30,60],state_window_indices=[len(s)-30,len(s)], registration_only=True,lag_identifiable=bool(case['scene']=='karman' and np.std(t[:,1])>1e-10)) data={} for role in ('target','zero','controlled'): with np.load(p/role/'fields.npz') as z:data[role]={k:z[k] for k in z.files} for role in ('zero','controlled'): name={'karman':'independent_phase8_diagnostic','steady':'late_snapshot','taylor':'common_event8'}[case['scene']] if case['scene']!='karman' or all(len(data[k]['ux'])==8 for k in (role,'target')): result['roles'][role][name]=compare_selected(data['target'],data[role]) else: result['roles'][role][name]=dict(status='unavailable_insufficient_periodic_crossings') result['roles'][role]['selection_metadata']=reports[role]['selection'] result['roles'][role]['field_window_tU_D']=reports[role]['field_window_tU_D'] if case['scene']=='karman' and all('cycle_mean_ux' in data[k] for k in (role,'target')): support=np.zeros(data[role]['cycle_mean_ux'].shape,bool);support[311:711,1360:2161]=True result['roles'][role]['E_complete_cycle_mean']=field_error( data[role]['cycle_mean_ux'],data[role]['cycle_mean_uy'], data['target']['cycle_mean_ux'],data['target']['cycle_mean_uy'],.005,support) if case['scene']=='taylor': clocks=[data[k]['absolute_steps'] for k in data] if not all(np.array_equal(clocks[0],v) for v in clocks[1:]):raise ValueError('event snapshot clocks differ') result['Taylor_injection']=dict(operator='all interior equilibrium rewrite; both buffers, identical event clock', strength=.000075,center=[600.,510.5],radius=80.,impulse_window='first two SI; retained native history', source_offsets=list(helpers.VORTEX_EVENT_OFFSETS)) result['same_clock_field_errors']=json.loads((p/'same_clock_errors.json').read_text()) run.write_json(p/'statistics.json',result) del data;gc.collect() print(json.dumps(dict(scene=case['scene'],statistics=str(p/'statistics.json')))) if __name__=='__main__':main()