131 lines
7.6 KiB
Python
131 lines
7.6 KiB
Python
"""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()
|