test: implement independently reviewed rising-crossing assessor v2

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,
+    )

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()

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)<len(cycles) and all(c['start']>=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())