#!/usr/bin/env python3 # SPDX-License-Identifier: MIT # Reproducible benchmark v1.5 - clock misalignment vs. real hydrogen loss on a dispenser line (illustrative model / 例示モデル / примерен модел) # 再現可能なベンチマーク v1.5：時刻ずれと実際の水素損失の区別（例示モデル） / Възпроизводим бенчмарк v1.5 (примерен модел) # Run: python3 mass_balance_clock.py > out.json (standard library only; output identical to mass_balance_clock.mjs) # Deterministic: own PRNG (mulberry32), own logarithm (series), only + - * / sqrt, identical evaluation order in JS. import json, math PARAMS = dict( version='1.5.0', seed=20261021, calSeed=20261022, bootSeed=20261023, dt=0.2, w=25, tail=3.0, # s, samples per window (W = w*dt = 5.0 s), s monitored after the fill nStations=50, fillsPerStation=8, # test: 400 fills per condition, instrument errors fixed per station calStations=100, calFillsPerStation=10, # calibration: 1000 fills per stratum, other stations, other seed targetFar=0.01, # nominal target: per-fill alarm probability on clean (no-leak) calibration fills calDeltas=[-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3], # declared clock-offset range |delta| <= 0.3 s: thresholds = worst case over these strata stage1Far=0.2, # B4 stage 1: B2 threshold at a 20 % per-fill rate (worst case over the strata) mpeMeter=0.015, # OIML R 139-1 Table 1, class 2, METER (per-meter value, not the system MPE) noiseAbs=0.15, noise='irwinhall4', arPhi=0.9, # g/s per sample (1 sd); noise model (assumption) qHold=[22, 39], rampUp=[5, 30], rampDown=[3, 10], hold=[60, 180], V=0.002, p0=[5, 20], pEnd=[70, 80], # m3 line volume between meters; MPa Tamb=[283.15, 303.15], Tpc=[233.15, 240.15], tauCool=[10, 40], tauWarm=200.0, # K, K, s, s pBias=0.005, TBias=1.0, tauT=5.0, tauMeter=0.3, # rel., K, s, s uTiming=[0.2, 0.4], # s, first-order timing allowances of B1 and B2 uLo=0.1, uHi=0.4, gq=1.0, gp=2.0, gm=2.0, # B3: u(t) = uLo + (uHi - uLo)*s(t); ramp scales g/s^2, MPa/s, g/s b5Shift=2, b5Rem=0.1, # B5: aligned shifts -2..+2 samples (+-0.4 s), first-order remainder for the half step (0.1 s) nb=5, sigM=0.1, kappa=0.5, rampGate=0.5, prior=[1.0, 0.05, 0.5, 0.3], # C/B4: 1 s blocks, misfit floor g/s, CUSUM allowance, steady-state gate, prior sd Rgas=4124.2, bAN=0.0076921, # J/(kg K); Abel-Noble co-volume m3/kg (Schefer et al. 2007) eos='abelNoble', # truth density model: 'abelNoble' or 'nistInterp' (model-mismatch test) lateProb=0.05, driftRate=0.003, tauOutLag=1.0, freezeMeter=10.0, freezeP=30.0, biasStep=0.02, jitter=0.1, mapLeaks=[0, 0.01, 0.02, 0.05], mapDeltas=[0, 0.2, 0.5, 1.0], mapFills=80, mapStride=5, timeRefErr=0.02, sensCal=500, sensFills=200, sensDeltas=[-0.3, 0, 0.3], bootB=1000, # NIST Chemistry WebBook (normal hydrogen, Leachman et al. 2009 EOS), density kg/m3 at T [K], p [MPa] nist=[[233.15, 5, 5.0294], [233.15, 20, 18.125], [233.15, 35, 28.621], [233.15, 70, 46.368], [233.15, 85, 52.104], [288.15, 5, 4.0833], [288.15, 20, 14.940], [288.15, 35, 23.995], [288.15, 70, 40.172], [288.15, 85, 45.630]], ) M32 = 0xFFFFFFFF LN2 = 0.6931471805599453 SQRT3 = math.sqrt(3) def imul(a, b): return (a * b) & M32 def mulberry32(seed): state = [seed & M32] def r(): a = (state[0] + 0x6D2B79F5) & M32 state[0] = a t = imul(a ^ (a >> 15), 1 | a) t = ((t + imul(t ^ (t >> 7), 61 | t)) & M32) ^ t return ((t ^ (t >> 14)) & M32) / 4294967296 return r def uni(r, lo, hi): return lo + (hi - lo) * r() def js_round(x): # JavaScript Math.round semantics return math.floor(x + 0.5) def ln(x): # natural log for x > 0 with only + - * / (identical in JS) e = 0 while x < 0.5: x = x * 2; e -= 1 while x >= 1: x = x / 2; e += 1 z = (x - 1) / (x + 1); z2 = z * z; s = 0.0; term = z for k in range(40): s = s + term / (2 * k + 1); term = term * z2 return 2 * s + e * LN2 def make_noise(r, model, phi): """returns a function giving the next standard-normal-like sample (mean 0, sd 1)""" if model == 'irwinhall4': # scaled sum of 4 uniforms (Irwin-Hall); light tails, |x| <= 2*sqrt(3) return lambda: (r() + r() + r() + r() - 2) * SQRT3 spare = [None] def polar(): # Marsaglia polar method: exact normal distribution if spare[0] is not None: v = spare[0]; spare[0] = None; return v while True: u = 2 * r() - 1; v = 2 * r() - 1; s = u * u + v * v if 0 < s < 1: break f = math.sqrt(-2 * ln(s) / s) spare[0] = v * f return u * f if model == 'gauss': return polar c = math.sqrt(1 - phi * phi); st = [None] def ar1(): # AR(1) with the same marginal sd: e[i] = phi*e[i-1] + sqrt(1-phi^2)*g g = polar() st[0] = g if st[0] is None else phi * st[0] + c * g return st[0] return ar1 def flow(t, f): # true outlet flow q(t) [g/s] if t <= 0: return 0 if t < f['tu']: return f['q'] * t / f['tu'] if t < f['tu'] + f['th']: return f['q'] d = t - f['tu'] - f['th'] if d < f['td']: return f['q'] * (1 - d / f['td']) return 0 def cum(t, f): # dispensed mass M(t) [g] q = f['q']; tu = f['tu']; th = f['th']; td = f['td'] if t <= 0: return 0 if t < tu: return q * t * t / (2 * tu) if t < tu + th: return q * tu / 2 + q * (t - tu) d = t - tu - th if d < td: return q * tu / 2 + q * th + q * (d - d * d / (2 * td)) return q * tu / 2 + q * th + q * td / 2 def rho(pMPa, T, P): # Abel-Noble real-gas density [kg/m3]: rho = p / (R T + b p) p = pMPa * 1e6 return p / (P['Rgas'] * T + P['bAN'] * p) def gen_station(r, P): m = P['mpeMeter'] return dict(bin=uni(r, -m, m), bout=uni(r, -m, m), bp=uni(r, -P['pBias'], P['pBias']), bT=uni(r, -P['TBias'], P['TBias'])) def gen_fill(r, P): return dict(q=uni(r, *P['qHold']), tu=uni(r, *P['rampUp']), td=uni(r, *P['rampDown']), th=uni(r, *P['hold']), p0=uni(r, *P['p0']), pEnd=uni(r, *P['pEnd']), Tamb=uni(r, *P['Tamb']), Tpc=uni(r, *P['Tpc']), tauC=uni(r, *P['tauCool']), uOn=r(), uFault=r(), uRef=r(), nseed=math.floor(r() * 4294967296)) def population(P, seed, nSt, perSt): r = mulberry32(seed); out = [] for s in range(nSt): st = gen_station(r, P) for _ in range(perSt): out.append((s, st, gen_fill(r, P))) return out def eos_corr(p, T, P): """relative density correction NIST/Abel-Noble - 1, linear in p (5 points) and T (2 isotherms), clamped; truth of the model-mismatch test only""" def c(q): return q[2] / rho(q[1], q[0], P) - 1 def at(r0): pts = P['nist'][r0:r0 + 5] if p <= pts[0][1]: return c(pts[0]) for j in range(4): a = pts[j]; b = pts[j + 1] if p <= b[1]: ca = c(a); cb = c(b) return ca + (cb - ca) * (p - a[1]) / (b[1] - a[1]) return c(pts[4]) c1 = at(0); c2 = at(5); T1 = P['nist'][0][0]; T2 = P['nist'][5][0] x = (T - T1) / (T2 - T1) if x < 0: x = 0 elif x > 1: x = 1 return c1 + (c2 - c1) * x def rho_true(p, T, P): r = rho(p, T, P) return r * (1 + eos_corr(p, T, P)) if P['eos'] == 'nistInterp' else r METHODS = ['A0', 'A', 'B1', 'B2', 'B3', 'B4', 'B5', 'C'] NS = 9 # record streams: 0..7 = METHODS (5 = B4, built after the fill), 8 = B5 without alignment (for "insufficient information") def solve_full(A, bvec): # Gauss-Jordan with partial pivoting: returns (x, inverse) or None k = len(bvec) M = [A[i][:] + [bvec[i]] + [1.0 if j == i else 0.0 for j in range(k)] for i in range(k)] for c in range(k): piv = c for rI in range(c + 1, k): if abs(M[rI][c]) > abs(M[piv][c]): piv = rI M[c], M[piv] = M[piv], M[c] pv = M[c][c] if abs(pv) < 1e-12: return None for j in range(2 * k + 1): M[c][j] = M[c][j] / pv for rI in range(k): if rI != c: fct = M[rI][c] for j in range(2 * k + 1): M[rI][j] = M[rI][j] - fct * M[c][j] return [M[i][k] for i in range(k)], [M[i][k + 1:] for i in range(k)] def gate(blk, c1): # records of the block stream among blocks whose B2 running maximum exceeds c1 mx = 0.0; out = [] for b in blk: if b[2] > c1 and b[0] > mx: mx = b[0]; out.append([b[0], b[1], b[3], b[4], b[5]]) return out def b4zmax(b4c, c1): for q in b4c: if q[0] > c1: return q[1] return 0.0 def zmax_of(rec): return rec[-1][0] if rec else 0.0 def first_alarm(rec, c): for r in rec: if r[0] > c: return r return None def simulate(st, f, cond, P, c1s=(), keep=False): """one fill; per method the records (new running maxima) [z, sign, sample, leak mass so far, undelivered mass]""" dt = P['dt']; w = P['w']; W = w * dt; nb = P['nb']; S = P['b5Shift'] n = js_round((f['tu'] + f['th'] + f['td'] + P['tail']) / dt) tflow = f['tu'] + f['th'] + f['td'] Mtot = cum(tflow + 1, f) aIn = dt / (P['tauMeter'] + dt); aOut = dt / (cond.get('tauOut', P['tauMeter']) + dt); aT = dt / (P['tauT'] + dt) aC = dt / (f['tauC'] + dt); aW = dt / (P['tauWarm'] + dt) noise = make_noise(mulberry32(f['nseed']), P['noise'], P['arPhi']) rl = mulberry32(f['nseed'] ^ 0x9E3779B9); rj = mulberry32(f['nseed'] ^ 0x85EBCA6B) delta = cond.get('delta', 0) drift = cond.get('drift', 0); late = cond.get('late', 0); leak = cond.get('leak', 0); fault = cond.get('fault', ''); jit = cond.get('jitter', 0) if P.get('timeRef') or cond.get('timeRef'): delta = delta - (delta + (2 * f['uRef'] - 1) * P['timeRefErr']) tOn = 0.0 if cond.get('leakFromStart') else f['uOn'] * (tflow - 5) tF = f['uFault'] * (tflow - 10) tFp = f['tu'] + f['uFault'] * (f['th'] - P['freezeP'] if f['th'] > P['freezeP'] else 0) sig = P['noiseAbs']; ub = P['mpeMeter'] / SQRT3; uB0 = P['uTiming'][0]; uB1 = P['uTiming'][1] nvar = sig * sig * 2 * w * dt * dt s2 = 2 * sig * sig / nb + P['sigM'] * P['sigM'] gqW = W * P['gq']; gpW = W * P['gp']; gmW = W * P['gm']; uLo = P['uLo']; uD = P['uHi'] - P['uLo']; rem = P['b5Rem'] kappa = P['kappa']; rgate = P['rampGate']; V = P['V']; eosN = P['eos'] == 'nistInterp' Tg = f['Tamb']; Ts = f['Tamb']; yIn = 0.0; yOut = 0.0; mPrev = 0.0; moPrev = 0.0; frozen = None; pFrozen = None Cin = [0.0] * (n + 1); Cout = [0.0] * (n + 1); mh = [0.0] * n; mo = [0.0] * n; pmA = [0.0] * n K = dict(mi=[0.0] * n, mo=mo, pm=pmA, mh=mh) if keep else None zmax = [0.0] * NS; rec = [[] for _ in range(NS)] by = 0.0; bq = 0.0; bm = 0.0; bc = 0; Qprev = None; Sp = 0.0; Sn = 0.0; ramp = 1.0 pr = P['prior'] Am = [[(1 / (pr[a] * pr[a]) if a == c else 0.0) for c in range(4)] for a in range(4)] Ainv = [[(pr[a] * pr[a] if a == c else 0.0) for c in range(4)] for a in range(4)] bv = [0.0, 0.0, 0.0, 0.0]; beta = [0.0, 0.0, 0.0, 0.0] blk = [] leakMass = 0.0 q = f['q']; pEnd = f['pEnd']; p0 = f['p0']; Tpc = f['Tpc']; Tamb = f['Tamb']; bout = st['bout']; bp = st['bp']; bT = st['bT'] for i in range(n): t = i * dt qo = flow(t, f) p = p0 + (pEnd - p0) * cum(t, f) / Mtot if i > 0: Tg = Tg + (aC * (Tpc - Tg) if qo > 0 else aW * (Tamb - Tg)) m = V * (rho_true(p, Tg, P) if eosN else rho(p, Tg, P)) * 1000 lk = leak * q * p / pEnd if (leak > 0 and t >= tOn and t < tflow) else 0 qi = qo + lk + ((m - mPrev) / dt if i > 0 else 0) mPrev = m d = delta + drift * t if jit > 0: d = d + jit * (2 * rj() - 1) # timestamp jitter: independent per sample, uniform +-jit xo = flow(t - d, f) yIn = yIn + aIn * (qi - yIn); yOut = yOut + aOut * (xo - yOut) bin_ = st['bin'] + (P['biasStep'] if (fault == 'biasStep' and t >= tOn) else 0) mi = yIn * (1 + bin_) + sig * noise() mout = yOut * (1 + bout) + sig * noise() if late > 0 and rl() < late: mout = moPrev # hold-last-value substitute for a late packet if fault == 'frozenMeter' and t >= tF and t < tF + P['freezeMeter']: if frozen is None: frozen = moPrev mout = frozen moPrev = mout pm = p * (1 + bp) if fault == 'frozenP' and t >= tFp and t < tFp + P['freezeP']: if pFrozen is None: pFrozen = pm pm = pFrozen Ts = Ts + aT * (Tg - Ts) mh[i] = V * rho(pm, Ts + bT, P) * 1000 # estimator always uses Abel-Noble mo[i] = mout; pmA[i] = pm Cin[i + 1] = Cin[i] + mi * dt; Cout[i + 1] = Cout[i] + mout * dt if keep: K['mi'][i] = mi if lk > 0: leakMass = leakMass + lk * dt undel = Mtot - cum(t, f) if i >= w: Sin = Cin[i + 1] - Cin[i + 1 - w]; Sout = Cout[i + 1] - Cout[i + 1 - w] U0 = math.sqrt((ub * Sin) * (ub * Sin) + (ub * Sout) * (ub * Sout) + nvar) dm = mh[i] - mh[i - w] R0 = Sin - Sout; R = R0 - dm dq = abs(mo[i] - mo[i - w]) s = max(dq / gqW, abs(pmA[i] - pmA[i - w]) / gpW, abs(dm) / gmW) if s > 1: s = 1 ramp = s u3 = uLo + uD * s; aR = abs(R) z = (abs(R0) / U0, aR / U0, (aR - uB0 * dq) / U0, (aR - uB1 * dq) / U0, (aR - u3 * dq) / U0) sg0 = 1 if R0 > 0 else -1; sg = 1 if R > 0 else -1 for k in range(5): if z[k] > zmax[k]: zmax[k] = z[k]; rec[k].append([z[k], sg0 if k == 0 else sg, i, leakMass, undel]) if i >= w + 2 * S: # B5: window ending S samples ago, outlet aligned by -S..S samples j = i - S Sin = Cin[j + 1] - Cin[j + 1 - w]; So0 = Cout[j + 1] - Cout[j + 1 - w] U = math.sqrt((ub * Sin) * (ub * Sin) + (ub * So0) * (ub * So0) + nvar) dmj = mh[j] - mh[j - w] zp = 0.0; zn = 0.0; zf = 0.0; r0 = 0.0 for s5 in range(-S, S + 1): So = Cout[j + 1 + s5] - Cout[j + 1 - w + s5] Rs = Sin - So - dmj dqs = abs(mo[j + s5] - mo[j - w + s5]) a = (Rs - rem * dqs) / U; b = (-Rs - rem * dqs) / U if s5 == -S or a < zp: zp = a if s5 == -S or b < zn: zn = b if s5 == 0: zf = (abs(Rs) - rem * dqs) / U; r0 = Rs z5 = zp if zp >= zn else zn; g5 = 1 if zp >= zn else -1 if z5 > zmax[6]: zmax[6] = z5; rec[6].append([z5, g5, i, leakMass, undel]) if zf > zmax[8]: zmax[8] = zf; rec[8].append([zf, 1 if r0 > 0 else -1, i, leakMass, undel]) if i >= 1: # C: 1 s blocks, recursive least squares on nuisance shapes, two one-sided CUSUMs md = (mh[i] - mh[i - 1]) / dt by = by + (mi - mout - md); bq = bq + mout; bm = bm + md; bc += 1 if bc == nb: Y = by / nb; Q = bq / nb; Mb = bm / nb D = 0.0 if Qprev is None else (Q - Qprev) / (nb * dt) Qprev = Q x = (1.0, Q, D, Mb) pred = 0.0; pv = 0.0 for a in range(4): pred = pred + x[a] * beta[a] for a in range(4): for c in range(4): pv = pv + x[a] * Ainv[a][c] * x[c] eps = (Y - pred) / math.sqrt(s2 + pv) if i >= w and ramp < rgate: Sp = max(0.0, Sp + eps - kappa); Sn = max(0.0, Sn - eps - kappa) for a in range(4): bv[a] = bv[a] + x[a] * Y / s2 for c in range(4): Am[a][c] = Am[a][c] + x[a] * x[c] / s2 sol = solve_full(Am, bv) if sol is not None: beta, Ainv = sol blk.append([Sp if Sp >= Sn else Sn, 1 if Sp >= Sn else -1, zmax[3], i, leakMass, undel]) by = 0.0; bq = 0.0; bm = 0.0; bc = 0 rec[7] = gate(blk, -1.0) rec4 = [gate(blk, c1) for c1 in c1s] rec[5] = rec4[0] if rec4 else [] suf = [0.0] * len(blk); mx = 0.0 for k in range(len(blk) - 1, -1, -1): if blk[k][0] > mx: mx = blk[k][0] suf[k] = mx b4c = [[blk[k][2], suf[k]] for k in range(len(blk)) if k == 0 or blk[k][2] > blk[k - 1][2]] return dict(rec=rec, rec4=rec4, b4c=b4c, tOn=tOn, leakMass=leakMass, n=n, keep=K, f=f, st=st) def wilson(x, n): z = 1.959963984540054; p = x / n; d = 1 + z * z / n c = (p + z * z / (2 * n)) / d; h = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d return [max(0, c - h), min(1, c + h)] def rr(v, k=4): s = 10 ** k return js_round(v * s) / s def solve(A, bvec): # Gauss-Jordan: returns (x, diag of inverse) k = len(bvec) M = [A[i][:] + [bvec[i]] + [1.0 if j == i else 0.0 for j in range(k)] for i in range(k)] for c in range(k): piv = c for rI in range(c + 1, k): if abs(M[rI][c]) > abs(M[piv][c]): piv = rI M[c], M[piv] = M[piv], M[c] pv = M[c][c] if abs(pv) < 1e-12: return None for j in range(2 * k + 1): M[c][j] = M[c][j] / pv for rI in range(k): if rI != c: fct = M[rI][c] for j in range(2 * k + 1): M[rI][j] = M[rI][j] - fct * M[c][j] return [M[i][k] for i in range(k)], [M[i][k + 1 + i] for i in range(k)] def classify(o, P, info, dref): """offline per-fill analysis: least squares of the residual on leak / meter-bias / clock-offset signatures""" K = o['keep']; n = o['n']; dt = P['dt'] mi = K['mi']; mo = K['mo']; pm = K['pm']; mh = K['mh'] qmax = 0.0 for i in range(n): if mo[i] > qmax: qmax = mo[i] nreg = 3 if info != 'I3' else 2 XtX = [[0.0] * nreg for _ in range(nreg)]; Xty = [0.0] * nreg; yy = 0.0; N = 0 w = P['w']; bx = [0.0] * nreg; by = 0.0; cnt = 0 # regression on non-overlapping block means (w samples) -> less autocorrelation for i in range(1, n - 1): der = (mo[i + 1] - mo[i - 1]) / (2 * dt) y = mi[i] - mo[i] xl = 1.0 if info != 'I1': y = y - (mh[i] - mh[i - 1]) / dt xl = pm[i] / 70 if info == 'I3': y = y - dref * der x = [xl, mo[i]] else: x = [xl, mo[i], der] for a in range(nreg): bx[a] = bx[a] + x[a] by = by + y; cnt += 1 if cnt == w: for a in range(nreg): bx[a] = bx[a] / w by = by / w for a in range(nreg): Xty[a] = Xty[a] + bx[a] * by for c in range(nreg): XtX[a][c] = XtX[a][c] + bx[a] * bx[c] yy = yy + by * by; N += 1 bx = [0.0] * nreg; by = 0.0; cnt = 0 sol = solve(XtX, Xty) if sol is None: return 'insufficient', False beta, inv = sol rss = yy for a in range(nreg): rss = rss - beta[a] * Xty[a] s2 = rss / (N - nreg) if rss > 0 else 0.0 seL = math.sqrt(s2 * inv[0]) if inv[0] > 0 else 0.0 Lmin = 0.01 * qmax offset = False if nreg == 3: seD = math.sqrt(s2 * inv[2]) if inv[2] > 0 else 0.0 offset = abs(beta[2]) > 1.96 * seD and abs(beta[2]) > 0.05 if 1.96 * seL > Lmin: return 'insufficient', offset return ('leak' if beta[0] > 1.96 * seL else 'noLeak'), offset def quantile_sorted(a, q): return a[min(len(a) - 1, math.floor(q * len(a)))] def qcut(srt, far): # alarm when z > cut: at most far*n of the fills return srt[len(srt) - math.floor(far * len(srt)) - 1] def fsum(a): s = 0 for x in a: s = s + x return s def boot_idx(P, nSt): r = mulberry32(P['bootSeed']) return [[math.floor(r() * nSt) for _ in range(nSt)] for _ in range(P['bootB'])] def boot_ci(BI, num, num2, den): """station-cluster bootstrap: resample whole stations (all 8 fills of a station together), percentile 95 % interval""" vals = [] for row in BI: a = 0; b = 0 for s in row: a += num[s] - (num2[s] if num2 else 0); b += den[s] vals.append(a / b) vals.sort(); B = len(vals) return [rr(vals[(B * 25) // 1000]), rr(vals[(B * 975) // 1000 - 1])] def stats(P, pop, outs, cond, thr, sset, nSt, BI): dt = P['dt']; leak = cond.get('leak', 0) > 0; out = {}; cnt = {} for k, mname in enumerate(METHODS): c = thr[k] alarms = neg = det = mis = miss = pre = ins = lowN = lowA = highN = highA = 0 stations = set(); ttd = []; lost = []; und = []; numA = [0] * nSt; numC = [0] * nSt; den = [0] * nSt for idx, o in enumerate(outs): sId = pop[idx][0]; den[sId] += 1 fa = first_alarm(o['rec4'][sset] if k == 5 else o['rec'][k], c); a = fa is not None if a: alarms += 1; stations.add(sId); numA[sId] += 1 if fa[1] < 0: neg += 1 und.append(fa[4]) if o['f']['q'] < 30.5: lowN += 1 if a: lowA += 1 else: highN += 1 if a: highA += 1 if leak: if not a: miss += 1; lost.append(o['leakMass']) elif fa[2] * dt < o['tOn']: pre += 1; lost.append(0.0) elif fa[1] > 0: det += 1; numC[sId] += 1; ttd.append(fa[2] * dt - o['tOn']); lost.append(fa[3]) else: mis += 1; lost.append(fa[3]) if k == 6 and not a and first_alarm(o['rec'][8], c) is not None: ins += 1 n = len(outs); lo, hi = wilson(alarms, n) e = dict(alarms=alarms, n=n, rate=rr(alarms / n), ci95=[rr(lo), rr(hi)]) if BI: e['ciBoot'] = boot_ci(BI, numA, None, den) e.update(negative=neg, stations=len(stations), nStations=nSt, lowFlow=[lowA, lowN], highFlow=[highA, highN], undeliveredMean=rr(fsum(und) / len(und), 1) if und else None) if leak: ttd.sort(); lost.sort() e.update(correct=det, mismatch=mis, preOnset=pre, missed=miss, correctRate=rr(det / n)) if BI: e['correctCiBoot'] = boot_ci(BI, numC, None, den) e.update(ttdMedian=rr(quantile_sorted(ttd, 0.5), 2) if ttd else None, ttdP90=rr(quantile_sorted(ttd, 0.9), 2) if ttd else None, lostMean=rr(fsum(lost) / len(lost), 2), lostP90=rr(quantile_sorted(lost, 0.9), 2)) if k == 6: e['insufficient'] = ins out[mname] = e; cnt[mname] = dict(numA=numA, numC=numC, den=den) return out, cnt def compact(ev): o = {} for m in METHODS: e = ev[m]; o[m] = dict(rate=e['rate']) if 'correct' in e: o[m].update(correct=e['correct'], preOnset=e['preOnset'], mismatch=e['mismatch'], missed=e['missed']) return o def calib(P, strata, deltas): """primary = worst case over the clock-offset strata (each stratum's (1-far) order statistic, then the maximum); cal0 = delta = 0 stratum only (as v1.4)""" i0 = deltas.index(0); far = P['targetFar'] def cut(arr, fr): return qcut(sorted(arr), fr) def wc(fz, fr): mx = None for sArr in strata: v = cut([fz(o) for o in sArr], fr) if mx is None or v > mx: mx = v return mx def zk(k): return lambda o: zmax_of(o['rec'][k]) c1 = wc(zk(3), P['stage1Far']); c10 = cut([zk(3)(o) for o in strata[i0]], P['stage1Far']) thr = [wc(lambda o: b4zmax(o['b4c'], c1), far) if k == 5 else wc(zk(k), far) for k in range(len(METHODS))] thr0 = [cut([b4zmax(o['b4c'], c10) for o in strata[i0]], far) if k == 5 else cut([zk(k)(o) for o in strata[i0]], far) for k in range(len(METHODS))] return dict(thr=thr, thr0=thr0, c1=c1, c10=c10) ROC_FAR = [0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5] def roc(P, calStrata, c1, testNeg, i0, sset, leakSets): """thresholds from calibration quantiles; FAR re-measured on separate negative TEST fills; leak fills: any alarm, first alarm after onset, correct diagnosis""" out = {} for k, m in enumerate(METHODS): def zc(o): return b4zmax(o['b4c'], c1) if k == 5 else zmax_of(o['rec'][k]) def zt(o): return zmax_of(o['rec4'][sset] if k == 5 else o['rec'][k]) cs = [sorted([zc(o) for o in sArr]) for sArr in calStrata]; ts = [[zt(o) for o in sArr] for sArr in testNeg] pts = [] for far in ROC_FAR: c = None for a in cs: v = qcut(a, far) if c is None or v > c: c = v fCal = 0 for a in cs: x = 0 for v in a: if v > c: x += 1 if x / len(a) > fCal: fCal = x / len(a) fWc = 0; f0 = 0 for si, a in enumerate(ts): x = 0 for v in a: if v > c: x += 1 r = x / len(a) if r > fWc: fWc = r if si == i0: f0 = r pt = dict(target=far, farCal=rr(fCal), farTest0=rr(f0), farTestWc=rr(fWc)) for lid, outs in leakSets: anyA = 0; after = 0; cor = 0 for o in outs: fa = first_alarm(o['rec4'][sset] if k == 5 else o['rec'][k], c) if fa is None: continue anyA += 1 if fa[2] * P['dt'] >= o['tOn']: after += 1 if fa[1] > 0: cor += 1 pt[lid] = [rr(anyA / len(outs)), rr(after / len(outs)), rr(cor / len(outs))] pts.append(pt) out[m] = pts return out MAP_ON = ['A', 'B2', 'B3', 'B4', 'B5', 'C'] def dmap(P, pop, thr, c1): """performance map of ONE specific offline classifier (whole-fill least squares, leak present from the start) plus the online detectors""" cells = []; IS = ('I1', 'I2', 'I3') sub = [pop[k * P['mapStride']] for k in range(P['mapFills'])] nSt = len(set(x[0] for x in sub)) for L in P['mapLeaks']: for dl in P['mapDeltas']: cnt = {I: dict(leak=0, noLeak=0, insufficient=0, offset=0) for I in IS} on = {m: [0, 0] for m in MAP_ON} for (_, st, f) in sub: o = simulate(st, f, dict(delta=dl, leak=L, leakFromStart=True), P, [c1], True) dref = dl + (2 * f['uRef'] - 1) * P['timeRefErr'] for I in IS: c, off = classify(o, P, I, dref) cnt[I][c] += 1; cnt[I]['offset'] += 1 if off else 0 for m in MAP_ON: k = METHODS.index(m); fa = first_alarm(o['rec'][k], thr[k]) if fa is not None: on[m][0] += 1 if fa[1] > 0: on[m][1] += 1 nn = len(sub) cell = dict(leak=L, delta=dl, n=nn, nStations=nSt) for I in IS: cell[I] = {k: rr(v / nn) for k, v in cnt[I].items()} cell['online'] = {m: [rr(on[m][0] / nn), rr(on[m][1] / nn)] for m in MAP_ON} cells.append(cell) return cells CONDITIONS = [ *[('d%s' % x, 'offset', dict(delta=x)) for x in [-1, -0.5, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.5, 1]], ('drift', 'timing', dict(drift=0.003)), ('jitter', 'timing', dict(jitter=0.1)), ('late', 'timing', dict(late=0.05)), ('lag', 'timing', dict(tauOut=1.0)), ('frozenMeter', 'fault', dict(fault='frozenMeter')), ('frozenP', 'fault', dict(fault='frozenP')), ('biasStep', 'fault', dict(fault='biasStep')), ('L1', 'leak', dict(leak=0.01)), *[('L%d_d%s' % (round(L * 100), x), 'leak', dict(leak=L, delta=x)) for L in [0.02, 0.05, 0.1] for x in [0, 0.2, 0.3, -0.3]], ] VARIANTS = [ ('base', 'sens', 'baseline at sensitivity sizes (500 cal fills per stratum, 200 test fills)', {}), ('gauss', 'sens', 'noise: true Gaussian (polar method)', dict(noise='gauss')), ('ar1', 'sens', 'noise: correlated AR(1), phi = 0.9', dict(noise='ar1')), ('w13', 'sens', 'window 2.6 s (13 samples)', dict(w=13)), ('w50', 'sens', 'window 10 s (50 samples)', dict(w=50)), ('mpe2', 'sens', 'per-meter error +-2 % (stress case; OIML R 139 gives 2 % as a SYSTEM MPE, not a meter value)', dict(mpeMeter=0.02)), ('mpe3', 'sens', 'per-meter error +-3 % (stress case; OIML R 139 gives 3 % as a SYSTEM in-service MPE, not a meter value)', dict(mpeMeter=0.03)), ('V5', 'sens', 'line volume 5 L', dict(V=0.005)), ('tauT15', 'sens', 'temperature-sensor lag 15 s', dict(tauT=15.0)), ('eos', 'mismatch', 'model mismatch: truth density interpolated from NIST (Leachman EOS), estimator Abel-Noble', dict(eos='nistInterp')), ('timeRef', 'upgrade', 'added measurement: independent time reference, residual error +-0.02 s', dict(timeRef=True)), ('tauT1', 'upgrade', 'added measurement: faster thermometer, lag 1 s', dict(tauT=1.0)), ('mpe05', 'upgrade', 'added measurement: more accurate flowmeters, +-0.5 % (hypothetical)', dict(mpeMeter=0.005)), ] SENS_CONDS = [('d0', dict(delta=0)), ('d0.2', dict(delta=0.2)), ('L2_d0', dict(leak=0.02)), ('L5_d0', dict(leak=0.05)), ('L2_d0.2', dict(leak=0.02, delta=0.2)), ('L5_d0.2', dict(leak=0.05, delta=0.2))] COMPARE = [ *[(cid, 'correct', a, b) for cid in ['L2_d0', 'L5_d0', 'L2_d0.2', 'L5_d0.2'] for a, b in [('B3', 'B2'), ('B4', 'B2'), ('B4', 'A'), ('B5', 'B2'), ('C', 'A'), ('C', 'B2')]], *[(cid, 'alarms', a, b) for cid in ['d0.2', 'd0.3', 'jitter'] for a, b in [('A', 'B2'), ('B3', 'B2'), ('B4', 'B2'), ('C', 'B2')]], ] COST_PI = [0.001, 0.01]; COST_R = [10, 100, 1000]; COST_NEG = ['d0', 'd0.2', 'd0.5', 'drift', 'lag'] def num(v): # print like JavaScript (1.0 -> 1) if isinstance(v, float) and v.is_integer(): return int(v) if isinstance(v, dict): return {k: num(x) for k, x in v.items()} if isinstance(v, (list, tuple)): return [num(x) for x in v] return v def run(P=PARAMS, progress=None): total = 1 + len(CONDITIONS) + len(VARIANTS) + 2; step = [0] def tick(): step[0] += 1 if progress: progress(step[0] / total) cal = population(P, P['calSeed'], P['calStations'], P['calFillsPerStation']) test = population(P, P['seed'], P['nStations'], P['fillsPerStation']) calStrata = [[simulate(st, f, dict(delta=d), P) for (_, st, f) in cal] for d in P['calDeltas']] cb = calib(P, calStrata, P['calDeltas']); tick() BI = boot_idx(P, P['nStations']) conds = []; keepOuts = {}; cnts = {} for cid, grp, c in CONDITIONS: outs = [simulate(st, f, c, P, [cb['c1'], cb['c10']]) for (_, st, f) in test] res, cnt = stats(P, test, outs, c, cb['thr'], 0, P['nStations'], BI) res0, _ = stats(P, test, outs, c, cb['thr0'], 1, P['nStations'], None) conds.append(dict(id=cid, group=grp, cond=c, res=res, res0=compact(res0))); cnts[cid] = cnt if grp == 'offset' and abs(c['delta']) <= 0.3 + 1e-9: keepOuts[cid] = outs if cid in ('L2_d0', 'L5_d0', 'L2_d0.2', 'L5_d0.2'): keepOuts[cid] = outs tick() comparisons = [] for cid, metric, a, b in COMPARE: ca = cnts[cid][a]; cbb = cnts[cid][b] na = ca['numA'] if metric == 'alarms' else ca['numC']; nb2 = cbb['numA'] if metric == 'alarms' else cbb['numC'] comparisons.append(dict(cond=cid, metric=metric, a=a, b=b, diff=rr((fsum(na) - fsum(nb2)) / fsum(ca['den'])), ciBoot=boot_ci(BI, na, nb2, ca['den']))) testNeg = [keepOuts['d%s' % d] for d in P['calDeltas']] leakSets = [(cid, keepOuts[cid]) for cid in ('L2_d0', 'L5_d0', 'L2_d0.2', 'L5_d0.2')] i0 = P['calDeltas'].index(0) rc = roc(P, calStrata, cb['c1'], testNeg, i0, 0, leakSets) rc0 = roc(P, [calStrata[i0]], cb['c10'], [keepOuts['d0']], 0, 1, leakSets) del calStrata, keepOuts sens = [] for vid, kind, name, ch in VARIANTS: Q = {**P, **ch} calQ = population(Q, Q['calSeed'], Q['calStations'], Q['calFillsPerStation'])[:Q['sensCal']] # regenerated with the variant's parameters, same seeds testQ = population(Q, Q['seed'], Q['nStations'], Q['fillsPerStation'])[:Q['sensFills']] strata = [[simulate(st, f, dict(delta=d), Q) for (_, st, f) in calQ] for d in Q['sensDeltas']] cq = calib(Q, strata, Q['sensDeltas']); nSt = Q['sensFills'] // Q['fillsPerStation'] rows = [] for cid, c in SENS_CONDS: c1s = [cq['c1'], cq['c10'], cb['c1']] if kind == 'mismatch' else [cq['c1'], cq['c10']] outs = [simulate(st, f, c, Q, c1s) for (_, st, f) in testQ] row = dict(id=cid, res=compact(stats(Q, testQ, outs, c, cq['thr'], 0, nSt, None)[0]), res0=compact(stats(Q, testQ, outs, c, cq['thr0'], 1, nSt, None)[0])) if kind == 'mismatch': row['nominal'] = compact(stats(Q, testQ, outs, c, cb['thr'], 2, nSt, None)[0]) rows.append(row) sens.append(dict(id=vid, kind=kind, name=name, thresholds=[rr(x) for x in cq['thr']], thresholds0=[rr(x) for x in cq['thr0']], rows=rows)); tick() mp = dmap(P, test, cb['thr'], cb['c1']); tick() CR = {c['id']: c['res'] for c in conds} cost = [] # expected cost per 1000 fills: (1 - pi) * P(unnecessary shutdown | negative condition) + pi * (r * P(not correctly diagnosed) + P(shutdown before onset)) for neg in COST_NEG: for L in [2, 5]: for pi in COST_PI: for r in COST_R: nres = CR[neg]; pos = CR['L%d_d0' % L] row = dict(neg=neg, leak=L / 100, pi=pi, ratio=r, cost={}) for m in METHODS: e = pos[m] row['cost'][m] = rr(1000 * ((1 - pi) * nres[m]['rate'] + pi * (((e['missed'] + e['mismatch']) / e['n']) * r + e['preOnset'] / e['n'])), 2) cost.append(row) dens = [dict(T=T, p=p, nist=d, abelNoble=rr(rho(p, T, P), 3), errPct=rr((rho(p, T, P) - d) / d * 100, 2)) for T, p, d in P['nist']] tick() return num({'benchmark': 'clock misalignment vs. real hydrogen loss (line-inventory mass balance)', 'version': P['version'], 'seed': P['seed'], 'calSeed': P['calSeed'], 'bootSeed': P['bootSeed'], 'methods': METHODS, 'calibration': dict(protocol='worst case over declared clock-offset strata', deltas=P['calDeltas'], nPerStratum=len(cal), targetFar=P['targetFar'], stage1Far=P['stage1Far'], thresholds=[rr(x) for x in cb['thr']], c1=rr(cb['c1']), thresholds0=[rr(x) for x in cb['thr0']], c10=rr(cb['c10'])), 'nTest': len(test), 'window': dict(samples=P['w'], seconds=rr(P['w'] * P['dt'], 3)), 'density': dens, 'conditions': conds, 'comparisons': comparisons, 'roc': rc, 'roc0': rc0, 'rocTargets': ROC_FAR, 'sensitivity': sens, 'map': mp, 'cost': cost}) if __name__ == '__main__': print(json.dumps(run(), indent=1, ensure_ascii=False))