diff --git a/tests/validation/cycle_v2/periodic_v2.py b/tests/validation/cycle_v2/periodic_v2.py new file mode 100644 index 0000000..474f169 --- /dev/null +++ b/tests/validation/cycle_v2/periodic_v2.py @@ -0,0 +1,102 @@ +"""Versioned ascending-zero DFG measurement; original data and gates unchanged.""" + +import numpy as np + +from periodic import GATES, REFERENCE_GATES, compare + +VERSION = "dfg-ascending-zero-v2" + + +def ascending_crossings(t, y): + t = np.asarray(t, dtype=float); y = np.asarray(y, dtype=float) + if t.ndim != 1 or y.shape != t.shape or not np.isfinite(t).all() or not np.isfinite(y).all(): + raise ValueError("Finite one-dimensional time/signal arrays required") + if np.any(np.diff(t) <= 0): + raise ValueError("Strictly increasing times required") + i = np.flatnonzero((y[:-1] < 0) & (y[1:] >= 0)) + return (t[i] - y[i] * (t[i+1] - t[i]) / (y[i+1] - y[i])).tolist() + + +def validate_rows(rows): + keys = ("tstar", "tstar_force", "cd", "cl", "dp", "mass", "flux_in", "flux_out") + values = np.asarray([[row[k] for k in keys] for row in rows], dtype=float) + if not np.isfinite(values).all(): + raise ValueError("Nonfinite observation") + if len(rows) and (np.any(np.diff(values[:, 0]) <= 0) or np.any(np.diff(values[:, 1]) <= 0)): + raise ValueError("Strictly increasing observation times required") + + +def extract(rows, min_time=250.0): + validate_rows(rows) + if len(rows) < 3: + return [] + t = np.array([r["tstar_force"] for r in rows]) + cl = np.array([r["cl"] for r in rows]) + ends = ascending_crossings(t, cl) + cycles = [] + for a, b in zip(ends, ends[1:]): + if a < min_time or not GATES["min_period"] < b - a < GATES["max_period"]: + continue + grid = np.linspace(a, b, 257) + c = dict(start=a, end=b, period=b - a, st=1 / (b - a)) + for key in ("cd", "cl", "dp", "mass", "flux_in", "flux_out"): + times = np.array( + [r["tstar" if key not in ("cd", "cl") else "tstar_force"] for r in rows] + ) + # No endpoint extrapolation: pressure/mass must cover the whole cycle. + if a < times[0] or b > times[-1]: + break + v = np.interp(grid, times, [r[key] for r in rows]) + c[key + "_min"] = float(v.min()) + c[key + "_max"] = float(v.max()) + c[key + "_midrange"] = float((v.min() + v.max()) / 2) + c[key + "_amplitude"] = float(v.max() - v.min()) + c[key + "_time_mean"] = float( + np.sum((v[1:] + v[:-1]) * np.diff(grid) / 2) / (b - a) + ) + if key in ("cd", "cl", "dp"): + c[key + "_waveform"] = v.tolist() + else: + c["flux_relative"] = abs( + c["flux_in_time_mean"] - c["flux_out_time_mean"] + ) / max(abs(c["flux_in_time_mean"]), 1e-30) + cycles.append(c) + return cycles + + +def assess(rows): + cycles = extract(rows, min_time=GATES["min_tstar"]) + tail = cycles[-GATES["cycles"] :] + checks = {} + checks["enough_cycles"] = len(tail) == GATES["cycles"] + if checks["enough_cycles"]: + checks["consecutive"] = all( + abs(a["end"] - b["start"]) < 1e-8 for a, b in zip(tail, tail[1:]) + ) + checks["recent"] = ( + 0 <= rows[-1]["tstar_force"] - tail[-1]["end"] <= GATES["max_period"] + ) + checks["shedding"] = ( + min(c["cl_amplitude"] for c in tail) >= GATES["min_lift_amplitude"] + ) + periods = [c["period"] for c in tail] + checks["period"] = ( + np.ptp(periods) / np.mean(periods) <= GATES["period_relative_range"] + ) + checks["extrema"] = all( + np.ptp([c[k + "_" + e] for c in tail]) <= GATES["extrema_absolute_range"] + for k in ("cd", "cl", "dp") + for e in ("min", "max") + ) + mass = [c["mass_time_mean"] for c in tail] + checks["mass"] = ( + np.ptp(mass) / abs(np.mean(mass)) <= GATES["mass_relative_range"] + ) + checks["flux"] = ( + max(c["flux_relative"] for c in tail) <= GATES["cycle_mean_flux_relative"] + ) + return dict( + periodic=bool(checks["enough_cycles"] and all(checks.values())), + checks={k: bool(v) for k, v in checks.items()}, + cycles=cycles, + ) diff --git a/tests/validation/cycle_v2/reanalyse.py b/tests/validation/cycle_v2/reanalyse.py new file mode 100644 index 0000000..94f258b --- /dev/null +++ b/tests/validation/cycle_v2/reanalyse.py @@ -0,0 +1,92 @@ +"""Versioned read-only DFG reanalysis; never overwrites a flow receipt.""" +import argparse +import csv +import hashlib +import io +import json +from pathlib import Path +import sys +import zipfile +import numpy as np +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +import periodic as frozen +import periodic_v2 as v2 + +OFFICIAL_SHA256 = '08db59605b2f60efc512d2bbc79628ea5d1cc720c0986938a16920884a0215be' +PREFIX = 'logs/results_cc2d_bench2_referencenew/' + + +def sha(path): return hashlib.sha256(Path(path).read_bytes()).hexdigest() + + +def official_reference(path): + if sha(path) != OFFICIAL_SHA256: raise ValueError('Wrong official archive') + with zipfile.ZipFile(path) as z: + force = np.loadtxt(io.BytesIO(z.read(PREFIX+'bdforces_q2_lv6_dt4'))) + pressure = np.loadtxt(io.BytesIO(z.read(PREFIX+'pointvalues_q2_lv6_dt4'))) + tf=force[:,1]*10; tp=pressure[:,1]*10 + if not np.isfinite(force).all() or not np.isfinite(pressure).all(): + raise ValueError('Nonfinite official input') + if np.any(np.diff(tf)<=0) or np.any(np.diff(tp)<=0): raise ValueError('Official time order') + # Entire supplied interval, restricted only to shared force/pressure coverage. + mask=(tf>=tp[0])&(tf<=tp[-1]); force=force[mask]; tf=tf[mask] + dp=np.interp(tf,tp,pressure[:,6]-pressure[:,11]) + rows=[dict(tstar=float(t),tstar_force=float(t),cd=float(f[3]),cl=float(f[4]),dp=float(p), + mass=1.,flux_in=1.,flux_out=1.) for t,f,p in zip(tf,force,dp)] + cycles=v2.extract(rows,min_time=0.) + if not cycles: raise ValueError('No complete official cycle') + ref=cycles[-1] + return ref,dict(archive_sha256=OFFICIAL_SHA256,selection='lv6/dt4 last complete ascending Cl=0 cycle in entire shared input coverage', + eligible_time=[float(tf[0]),float(tf[-1])],complete_cycles=len(cycles), + pressure='columns7 minus12, interpolated at original force times; physical time multiplied by10', + artificial_fields='mass/flux placeholders only; official archive is not used to certify conservation') + + +def main(): + parser=argparse.ArgumentParser(description=__doc__) + parser.add_argument('--official-archive',type=Path,required=True) + parser.add_argument('--history',action='append',required=True,help='Unique label=CSV path; all cases reported') + parser.add_argument('--output',type=Path,required=True) + a=parser.parse_args() + if a.output.exists(): raise ValueError('Use a new output directory') + cases=[]; labels=set() + reference,reference_info=official_reference(a.official_archive) + old_ref_path=Path(frozen.__file__).parent/'data/dfg_re100/reference_cycle.json' + old_ref=json.loads(old_ref_path.read_text()) + for item in a.history: + label,path=item.split('=',1); path=Path(path) + if label in labels: raise ValueError('Duplicate case label') + labels.add(label) + rows=[{k:float(v) for k,v in r.items()} for r in csv.DictReader(path.open())] + source_receipt_path=path.with_name('summary.json') + source_receipt=json.loads(source_receipt_path.read_text()) + if rows != source_receipt.get('history'): + raise ValueError('CSV differs from original summary history') + old=frozen.assess(rows); new=v2.assess(rows) + for key in ('periodic', 'checks', 'cycles'): + if old[key] != source_receipt.get(key): + raise ValueError('Frozen reassessment differs from original summary') + old_comparison=frozen.compare(old['cycles'][-1],old_ref) if old['cycles'] else None + new_comparison=v2.compare(new['cycles'][-1],reference) if new['cycles'] else None + # This is assessment of recorded observables, not proof of unrecorded population health. + healthy=bool(rows and all(.5<=r['rho_min']<=r['rho_max']<=1.5 and 0<=r['max_mach']<=.2 for r in rows)) + cases.append(dict(label=label,input_csv_sha256=sha(path), + original_summary_sha256=sha(source_receipt_path), + original_identity={k:source_receipt[k] for k in ('case','collision','observation_mode','numerical_source_commit','checkout','checkout_revision','mode','status') if k in source_receipt}, + observations=len(rows), + recorded_end_tstar=rows[-1]['tstar'] if rows else None,recorded_observations_healthy=healthy, + frozen=old,frozen_reference=old_comparison,v2=new,v2_reference=new_comparison, + v2_recorded_observable_agreement=bool(healthy and new['periodic'] and new_comparison and new_comparison['reference_agreement']))) + a.output.mkdir(parents=True) + result=dict(assessment_version=v2.VERSION,original_receipts_unchanged=True, + phase='Every ascending Cl<0 to Cl>=0 crossing; linear interpolation, no suppression', + gates=v2.GATES,reference_gates=v2.REFERENCE_GATES,reference=reference,reference_provenance=reference_info, + code_sha256={Path(p).name:sha(p) for p in [__file__,v2.__file__,frozen.__file__]}, + frozen_reference_sha256=sha(old_ref_path),cases=cases) + (a.output/'reanalysis.json').write_text(json.dumps(result,indent=2,allow_nan=False)+'\n') + print(json.dumps([dict(case=c['label'],old_periodic=c['frozen']['periodic'], + v2_checks=c['v2']['checks'],v2_reference=c['v2_reference'], + recorded_observable_agreement=c['v2_recorded_observable_agreement']) for c in cases],indent=2)) + + +if __name__=='__main__': main() diff --git a/tests/validation/cycle_v2/test_periodic_v2.py b/tests/validation/cycle_v2/test_periodic_v2.py new file mode 100644 index 0000000..a5aed8f --- /dev/null +++ b/tests/validation/cycle_v2/test_periodic_v2.py @@ -0,0 +1,88 @@ +"""Predeclared synthetic measurement tests; no flow outcome is used for tuning.""" +from pathlib import Path +import sys +import unittest +import numpy as np +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +import periodic +import periodic_v2 as v2 + + +def signal(periods=None, noise=0., offset=0., amplitude=1., irregular=False, harmonic=0.): + periods = np.full(20, 3.3) if periods is None else np.asarray(periods) + boundaries = 250.37 + np.r_[0., np.cumsum(periods)] + t = np.arange(boundaries[0]-.1, boundaries[-1]+.1, .02) + if irregular: t = t[np.arange(len(t)) % 7 != 0] + phase = np.interp(t, boundaries, np.arange(len(boundaries))) * 2*np.pi + lift = amplitude*(np.sin(phase) + harmonic*np.sin(2*phase)) + offset + # Fixed seed and fixed peak-local mask; not fitted to any flow dataset. + lift += noise*np.random.default_rng(1729).uniform(-1,1,len(t))*(abs(np.sin(phase))>.9) + return [dict(tstar=x+.01, tstar_force=x, cl=l, + cd=3.2+.03*np.cos(2*p), dp=2.45+.05*np.cos(p), + mass=1., flux_in=1., flux_out=1.) for x,l,p in zip(t,lift,phase)] + + +def test_constant_period_peak_noise_and_mean_offset(): + for noise, offset in [(0.,0.),(.0001,0.),(.001,0.),(.001,.08)]: + result=v2.assess(signal(noise=noise,offset=offset)) + assert result['periodic'], result['checks'] + periods=np.array([c['period'] for c in result['cycles'][-5:]]) + assert max(abs(periods-3.3))<1e-5 + + +def test_noncommensurate_asymmetric_waveform(): + result=v2.assess(signal(periods=[3.313]*20,harmonic=.1,offset=.08,noise=.001)) + assert result['periodic'], result['checks'] + periods=np.array([c['period'] for c in result['cycles'][-5:]]) + assert max(abs(periods-3.313))<5e-5 + + +def test_true_drift_is_rejected(): + result=v2.assess(signal(periods=np.arange(20)*.01+3.2)) + assert not result['periodic'] and not result['checks']['period'] + measured=np.array([c['period'] for c in result['cycles'][-5:]]) + assert np.ptp(measured)/np.mean(measured)>.005 + + +def test_insufficient_cycles_and_amplitude(): + assert not v2.assess(signal(periods=[3.3]*4))['periodic'] + result=v2.assess(signal(amplitude=.05)) + assert not result['periodic'] and not result['checks']['shedding'] + + +def test_exact_zero_plateau_and_crossing_noise(): + assert v2.ascending_crossings([0,1,2,3,4,5],[-1,0,0,1,-1,0])==[1.,5.] + # All crossings retained, even a short sign reversal: no hidden debounce. + assert len(v2.ascending_crossings([0,1,2,3,4,5],[-1,.001,-.001,.001,1,2]))==2 + assert v2.ascending_crossings([0,1,2],[1,.1,2])==[] + + +def test_extra_and_missing_cycles_not_reconnected(): + for last in [[3.3,.2,3.1,3.3],[3.3,6.6,3.3,3.3]]: + result=v2.assess(signal(periods=[3.3]*12+last)) + assert not result['periodic'] and not result['checks']['consecutive'] + + +def test_invalid_inputs_and_irregular_sampling(): + for t,y in [([0,0],[-1,1]),([1,0],[-1,1]),([0,1],[-1,float('nan')])]: + with unittest.TestCase().assertRaises(ValueError):v2.ascending_crossings(t,y) + rows=signal();rows[1]['dp']=float('inf') + with unittest.TestCase().assertRaises(ValueError):v2.assess(rows) + assert v2.assess(signal(irregular=True))['periodic'] + + +def test_coverage_and_unchanged_gates(): + rows=signal() + cycles=v2.extract(rows,min_time=250.) + # Move end-sample coverage late; force remains unchanged. Partial cycles excluded. + delayed=[dict(r,tstar=r['tstar']+10.) for r in rows] + other=v2.extract(delayed,min_time=250.) + assert len(other)=delayed[0]['tstar'] for c in other) + assert v2.GATES is periodic.GATES and v2.REFERENCE_GATES is periodic.REFERENCE_GATES + assert v2.compare is periodic.compare + + +if __name__=='__main__': + result=unittest.TextTestRunner(verbosity=2).run(unittest.TestSuite( + unittest.FunctionTestCase(fn) for name,fn in sorted(globals().items()) if name.startswith('test_'))) + raise SystemExit(not result.wasSuccessful())