feat: implement Pillars 2-4 and fix a production false-positive in Pillar 1

Pillar 1 v1.1 (shipped earlier) flagged benign browsing bursts as
beaconing: on short captures, grid periods longer than the observed
span never complete a cycle, so significance saturates to 1 and the
score sticks at exactly 0.5. v1.2 caps the grid at a third of the
observed span, requires a minimum number of completed cycles, and
excludes protocols that are periodic by nature (NTP, DNS, mDNS,
ICMP). Validated against 300 seeds per scenario: false positives on
bursty traffic drop from the 290-300/300 that v1.1 produced to 0/300,
with every beacon-detection scenario scoring unchanged.

Add Pillar 2 (persistent/long-lived sessions, with Zeek's fragmented
TCP sessions rejoined after a 5-minute gap), Pillar 3 (DNS tunneling
via Shannon-entropy information volume per subdomain per hour; the
Jensen-Shannon divergence was measured and does not separate tunnels
from benign CDN subdomains, so it stays in the evidence at zero
weight instead of in the score), and Pillar 4 (deterministic
noisy-OR fusion across the three pillars).

The fusion formula was validated against 8 labeled real-traffic
vectors (IoT-23, CTU-Normal) with tolerance <= 1e-4. A known
limitation is kept, not hidden: a benign Philips HUE bulb's strict
periodic polling and a 21.9h session score it above three real C2
hosts. No weight or threshold was tuned to mask this; it is
documented in the architecture docs instead.

Also fixes real bugs found along the way: ingest_dir's partition
count summed file names instead of distinct dates, DuckDB errors
were counted as bad rows, /api/alerts always returned an empty list
regardless of the last analysis, and the dashboard still claimed
every pillar was unimplemented.

Restores four design corrections that were missing from
docs/architecture.md (KL to Jensen-Shannon, the persistent-session
detector, the JA4 licensing boundary, and deterministic reduction
order), and catalogs the public datasets used for validation
(IoT-23, CTU-Normal, Mendeley DTQBC) with their licenses and SHA-256
checksums under research/validation/.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VAyr8d7HuAdUGB6aA5xmDD

research/validation/datasets.tsv

@@ -0,0 +1,10 @@
+dataset	url	license	file	sha256
+IoT-23 34-1 (Mirai)	https://mcfp.felk.cvut.cz/publicDatasets/IoT-23-Dataset/IndividualScenarios/CTU-IoT-Malware-Capture-34-1/	CC BY 2.0	2018-12-21-15-50-14-192.168.1.195.pcap	92ec7e2f6658ee4b007d0b816986c46cc0338bc5e2bec6ceaaca566c695e4699
+IoT-23 20-1 (Torii)	https://mcfp.felk.cvut.cz/publicDatasets/IoT-23-Dataset/IndividualScenarios/CTU-IoT-Malware-Capture-20-1/	CC BY 2.0	2018-10-02-13-12-30-192.168.100.103.pcap	4d0a00c33a4bf11228158baddd47c38f702731ea8d75ef1201061fc8bbd878b9
+IoT-23 8-1 (Hakai)	https://mcfp.felk.cvut.cz/publicDatasets/IoT-23-Dataset/IndividualScenarios/CTU-IoT-Malware-Capture-8-1/	CC BY 2.0	2018-07-31-15-15-09-192.168.100.113.pcap	80dcc2602519479ddcde889fa902fee19a76696630811452f8df38888af894f2
+IoT-23 42-1 (Trojan)	https://mcfp.felk.cvut.cz/publicDatasets/IoT-23-Dataset/IndividualScenarios/CTU-IoT-Malware-Capture-42-1/	CC BY 2.0	2019-01-10-14-34-38-192.168.1.197.pcap	7573797f17aae96e3804a3cfb22d531b74c17bbae740de1e14048f89130fba37
+IoT-23 4-1 (Philips HUE)	https://mcfp.felk.cvut.cz/publicDatasets/IoT-23-Dataset/IndividualScenarios/CTU-Honeypot-Capture-4-1/	CC BY 2.0	2018-10-25-14-06-32-192.168.1.132.pcap.xz	bc5a57306d9a5f614b09d9c0976f007667fee5b6de3427fbab094fa3084b9b64
+CTU-Normal-20	https://mcfp.felk.cvut.cz/publicDatasets/CTU-Normal-20/	unspecified	2017-04-30_win-normal.pcap	c78a00728f115ab8fc495fd395894c43fd90cb61054eb9ef4f2d024658d5ffc5
+CTU-Normal-18	https://mcfp.felk.cvut.cz/publicDatasets/CTU-Normal-18/	unspecified	2017-09-05_capture-win4.pcap	e09df1ac62d9753ed2360373915427432292192dc06e20572e45716a2932c8ce
+CTU-Normal-26	https://mcfp.felk.cvut.cz/publicDatasets/CTU-Normal-26/	unspecified	2017-04-28_normal.pcap	b66103fd89c8f0fc53b76e718f53320f3af27edd302c58695fc1870ede3663c5
+Mendeley dtqbc-m-test	https://data.mendeley.com/public-api/datasets/mzn9hvdcxg	CC BY 4.0	dtqbc-m-test.csv	48f109f501e01696e78d22623e4f7938784071219acbf8924ea2211f3473a4f8

research/validation/mkpcap_idle.py

@@ -0,0 +1,35 @@
+# Writes a pcap: one TCP session 10.0.0.5:40000 -> 203.0.113.9:443 with a 600 s idle gap
+# (longer than Zeek's tcp_inactivity_timeout = 5 min), then a fresh connection from a new port.
+import struct, sys
+out = open(sys.argv[1], 'wb')
+out.write(struct.pack('<IHHiIII', 0xa1b2c3d4, 2, 4, 0, 0, 65535, 1))
+def ip(a): return bytes(int(x) for x in a.split('.'))
+def pkt(t, src, sp, dst, dp, flags, seq, ack, payload=b''):
+    tcp = struct.pack('!HHIIBBHHH', sp, dp, seq, ack, 5 << 4, flags, 65535, 0, 0) + payload
+    iph = struct.pack('!BBHHHBBH4s4s', 0x45, 0, 20 + len(tcp), 1, 0, 64, 6, 0, ip(src), ip(dst))
+    frame = b'\x00\x11\x22\x33\x44\x55' + b'\x66\x77\x88\x99\xaa\xbb' + b'\x08\x00' + iph + tcp
+    sec = int(t); usec = int(round((t - sec) * 1e6))
+    out.write(struct.pack('<IIII', sec, usec, len(frame), len(frame)) + frame)
+SYN, ACK, PSH, FIN = 0x02, 0x10, 0x08, 0x01
+C, S = '10.0.0.5', '203.0.113.9'
+t0 = 1793400000.0
+cs, ss = 1000, 5000
+pkt(t0, C, 40000, S, 443, SYN, cs, 0); cs += 1
+pkt(t0+0.05, S, 443, C, 40000, SYN|ACK, ss, cs); ss += 1
+pkt(t0+0.10, C, 40000, S, 443, ACK, cs, ss)
+# heartbeat every 60 s for 30 min (stays alive)
+t = t0 + 1
+for i in range(30):
+    pkt(t, C, 40000, S, 443, PSH|ACK, cs, ss, b'x'*20); cs += 20
+    pkt(t+0.05, S, 443, C, 40000, PSH|ACK, ss, cs, b'y'*40); ss += 40
+    t += 60
+# idle 600 s, then resume heartbeats for 30 min more, then FIN
+t += 600
+for i in range(30):
+    pkt(t, C, 40000, S, 443, PSH|ACK, cs, ss, b'x'*20); cs += 20
+    pkt(t+0.05, S, 443, C, 40000, PSH|ACK, ss, cs, b'y'*40); ss += 40
+    t += 60
+pkt(t, C, 40000, S, 443, FIN|ACK, cs, ss); cs += 1
+pkt(t+0.05, S, 443, C, 40000, FIN|ACK, ss, cs); ss += 1
+pkt(t+0.10, C, 40000, S, 443, ACK, cs, ss)
+out.close()

research/validation/p1real.py

@@ -0,0 +1,45 @@
+# Applies Pillar 1 v1.1 and v1.2 (and ablations) to real Zeek conn.log files; tuple label from IoT-23 labeled file.
+import sys, os, collections, importlib.util
+here = os.path.dirname(os.path.abspath(__file__))
+spec = importlib.util.spec_from_file_location('sim', os.path.join(here, 'p1v12sim_lib.py')); sim = importlib.util.module_from_spec(spec); spec.loader.exec_module(sim)
+def read(path):
+    fields = None
+    for line in open(path, errors='replace'):
+        line = line.rstrip('\n')
+        if line.startswith('#fields'): fields = line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        yield dict(zip(fields, line.split('\t')))
+def labels(path):
+    lab = {}
+    if not path or not os.path.exists(path): return lab
+    for r in read(path):
+        last = ' '.join(v for k, v in r.items() if k in ('tunnel_parents', 'label', 'detailed-label', 'det_label'))
+        mal = 'Malicious' in last or 'Malicious' in str(list(r.values())[-3:])
+        k = (r['id.orig_h'], r['id.resp_h'], r['id.resp_p'])
+        lab[k] = lab.get(k, False) or mal
+    return lab
+import ipaddress
+EXCL = {('udp','5355'),('udp','546'),('udp','547'),('udp', '123'), ('udp', '53'), ('tcp', '53'), ('udp', '5353'), ('udp', '137'), ('udp', '138'), ('udp', '67'), ('udp', '68'), ('udp', '1900')}
+conn, labf = sys.argv[1], (sys.argv[2] if len(sys.argv) > 2 else None)
+lab = labels(labf)
+tup = collections.defaultdict(list); protos = {}
+for r in read(conn):
+    k = (r['id.orig_h'], r['id.resp_h'], r['id.resp_p']); tup[k].append(float(r['ts'])); protos.setdefault(k, set()).add(r['proto'])
+rows = []
+for k, ts in sorted(tup.items()):
+    if len(ts) < 20: continue
+    a = sim.v11(ts)[0]; b = sim.v12(ts)
+    span_only = sim.v12(ts, gate=1e-12)[0]; gate_only = sim.v12(ts, cycles=1e-12)[0]
+    try:
+        a_ = ipaddress.ip_address(k[1]); mc = a_.is_multicast or k[1] == '255.255.255.255'
+    except ValueError:
+        mc = False
+    excl = mc or any(p not in ('tcp', 'udp') for p in protos[k]) or any((p, k[2]) in EXCL for p in protos[k])
+    rows.append((k, len(ts), a, span_only, gate_only, b[0], b[1], b[2], excl, lab.get(k)))
+print(f"== {conn}: tuples with n>=20: {len(rows)}")
+print("tuple n v1.1 span_only gate_only v1.2 status period excluded label")
+for r in rows:
+    if max(r[2], r[5]) >= .40:
+        print(f"  {r[0]} n={r[1]} v11={r[2]:.3f} span={r[3]:.3f} gate={r[4]:.3f} v12={r[5]:.3f} {r[6]} per={r[7] and round(r[7])} excl={r[8]} mal={r[9]}")
+c = lambda i, ex: sum(1 for r in rows if r[i] >= .4 and (not ex or not r[8]))
+print(f"  alerts>=0.40: v1.1={c(2,False)} span_only={c(3,False)} gate_only={c(4,False)} v1.2={c(5,False)} v1.2+exclusion={c(5,True)}")

research/validation/p1v12sim.py

@@ -0,0 +1,104 @@
+# Pillar 1: v1.1 (implemented) vs v1.2 (proposed) on synthetic scenarios, 300 seeds each.
+import math, random, statistics, sys
+GRID = [600 * 72 ** (j / 199) for j in range(200)]
+TWO_PI = 2 * math.pi
+
+def qcont(xs, q):
+    xs = sorted(xs); h = (len(xs) - 1) * q; lo = math.floor(h); hi = math.ceil(h)
+    return xs[lo] + (h - lo) * (xs[hi] - xs[lo])
+
+def core(ts, grid):
+    t0 = ts[0]; n = len(ts); Ps = []
+    for p in grid:
+        f = 1 / p; c = s = 0.0
+        for t in ts:
+            ph = TWO_PI * f * (t - t0); c += math.cos(ph); s += math.sin(ph)
+        Ps.append((c * c + s * s) / n)
+    return Ps
+
+def disp(ts):
+    d = [b - a for a, b in zip(ts, ts[1:])]
+    med = qcont(d, .5)
+    if med <= 1e-9: return None, med
+    madm = qcont([abs(x - med) for x in d], .5)
+    return madm / med, med
+
+def v11(ts):
+    ts = sorted(ts)
+    if len(ts) < 20: return 0.0, 'insufficient_data', None
+    Ps = core(ts, GRID); pmax = max(Ps)
+    sig = min(1, max(0, (pmax - math.log(200)) / (6 * math.log(10))))
+    dn, med = disp(ts)
+    if dn is None: return 0.5 * sig, 'degenerate_deltas', None
+    reg = min(1, max(0, 1 - 2 * dn))
+    per = max(p for p, v in zip(GRID, Ps) if v >= .5 * pmax)
+    return 0.5 * sig + 0.5 * reg, 'ok', per
+
+def v12(ts, cycles=3, gate=0.25):
+    ts = sorted(ts)
+    if len(ts) < 20: return 0.0, 'insufficient_data', None
+    span = ts[-1] - ts[0]
+    grid = [p for p in GRID if p <= span / cycles]
+    if not grid: return 0.0, 'insufficient_span', None
+    Ps = core(ts, grid); pmax = max(Ps); M = len(grid)
+    sig = min(1, max(0, (pmax - math.log(M)) / (6 * math.log(10))))
+    dn, med = disp(ts)
+    if dn is None: return 0.0, 'degenerate_deltas', None
+    reg = min(1, max(0, 1 - 2 * dn))
+    g = min(1.0, reg / gate)
+    per = max(p for p, v in zip(grid, Ps) if v >= .5 * pmax)
+    return 0.5 * sig * g + 0.5 * reg, 'ok', per
+
+B = 1_790_000_000.0; T0 = 14400
+def perfect(r): return [B + k * T0 for k in range(42)]
+def grid(r): return [B + k * T0 + T0 * r.uniform(-.18, .18) for k in range(42)]
+def walk(r):
+    t, o = B, []
+    for _ in range(42): o.append(t); t += T0 * (1 + r.uniform(-.18, .18))
+    return o
+def poisson(r):
+    t, o = B, []
+    while len(o) < 40: o.append(t); t += r.expovariate(40 / (7 * 86400))
+    return o
+def dropout(r):  # 4 h beacon, +-10% cumulative jitter, 30% of beats missed (host asleep)
+    t, o = B, []
+    for _ in range(60):
+        if r.random() > .3: o.append(t)
+        t += T0 * (1 + r.uniform(-.1, .1))
+    return o
+def bursty(r, span):  # browsing-like: Poisson cluster starts, geometric cluster size, ~2 s intra-gaps
+    o = []; t = B
+    while True:
+        t += r.expovariate(1 / 900)
+        if t > B + span: break
+        k = 1
+        while r.random() < .8: k += 1
+        u = t
+        for _ in range(k): o.append(u); u += r.expovariate(1 / 2)
+    return o
+def bursty3h(r): return bursty(r, 3 * 3600)
+def bursty24h(r): return bursty(r, 86400)
+def fast_short(r):  # 15 min beacon observed only 24 h, 10% cumulative jitter
+    t, o = B, []
+    while t < B + 86400: o.append(t); t += 900 * (1 + r.uniform(-.1, .1))
+    return o
+
+def q(xs):
+    xs = sorted(xs); n = len(xs)
+    return f"min={xs[0]:.3f} p05={xs[n//20]:.3f} med={statistics.median(xs):.3f} p95={xs[-max(1,n//20)]:.3f} max={xs[-1]:.3f}"
+seeds = range(1, int(sys.argv[1]) + 1 if len(sys.argv) > 1 else 301)
+for name, gen in [('perfect', perfect), ('grid', grid), ('walk', walk), ('dropout', dropout), ('fast_short', fast_short),
+                  ('poisson', poisson), ('bursty3h', bursty3h), ('bursty24h', bursty24h)]:
+    a, b, pa, pb, st = [], [], [], [], {}
+    for s in (seeds if name != 'perfect' else [1]):
+        ts = gen(random.Random(s))
+        x = v11(ts); y = v12(ts)
+        a.append(x[0]); b.append(y[0]); st[y[1]] = st.get(y[1], 0) + 1
+        if x[2]: pa.append(x[2])
+        if y[2]: pb.append(y[2])
+    print(f"{name:10s} v1.1 {q(a)} | >=0.40: {sum(v >= .4 for v in a)}/{len(a)}")
+    print(f"{name:10s} v1.2 {q(b)} | >=0.40: {sum(v >= .4 for v in b)}/{len(b)} status={st}")
+    if name in ('grid', 'walk', 'dropout'):
+        print(f"{'':10s} period within 10% of {T0}: v1.1 {sum(abs(p-T0)<=.1*T0 for p in pa)}/{len(pa)} v1.2 {sum(abs(p-T0)<=.1*T0 for p in pb)}/{len(pb)}")
+    if name == 'fast_short':
+        print(f"{'':10s} period within 10% of 900: v1.2 {sum(abs(p-900)<=90 for p in pb)}/{len(pb)}")

research/validation/p1v12sim_lib.py

@@ -0,0 +1,51 @@
+# Pillar 1: v1.1 (implemented) vs v1.2 (proposed) on synthetic scenarios, 300 seeds each.
+import math, random, statistics, sys
+GRID = [600 * 72 ** (j / 199) for j in range(200)]
+TWO_PI = 2 * math.pi
+
+def qcont(xs, q):
+    xs = sorted(xs); h = (len(xs) - 1) * q; lo = math.floor(h); hi = math.ceil(h)
+    return xs[lo] + (h - lo) * (xs[hi] - xs[lo])
+
+def core(ts, grid):
+    t0 = ts[0]; n = len(ts); Ps = []
+    for p in grid:
+        f = 1 / p; c = s = 0.0
+        for t in ts:
+            ph = TWO_PI * f * (t - t0); c += math.cos(ph); s += math.sin(ph)
+        Ps.append((c * c + s * s) / n)
+    return Ps
+
+def disp(ts):
+    d = [b - a for a, b in zip(ts, ts[1:])]
+    med = qcont(d, .5)
+    if med <= 1e-9: return None, med
+    madm = qcont([abs(x - med) for x in d], .5)
+    return madm / med, med
+
+def v11(ts):
+    ts = sorted(ts)
+    if len(ts) < 20: return 0.0, 'insufficient_data', None
+    Ps = core(ts, GRID); pmax = max(Ps)
+    sig = min(1, max(0, (pmax - math.log(200)) / (6 * math.log(10))))
+    dn, med = disp(ts)
+    if dn is None: return 0.5 * sig, 'degenerate_deltas', None
+    reg = min(1, max(0, 1 - 2 * dn))
+    per = max(p for p, v in zip(GRID, Ps) if v >= .5 * pmax)
+    return 0.5 * sig + 0.5 * reg, 'ok', per
+
+def v12(ts, cycles=3, gate=0.25):
+    ts = sorted(ts)
+    if len(ts) < 20: return 0.0, 'insufficient_data', None
+    span = ts[-1] - ts[0]
+    grid = [p for p in GRID if p <= span / cycles]
+    if not grid: return 0.0, 'insufficient_span', None
+    Ps = core(ts, grid); pmax = max(Ps); M = len(grid)
+    sig = min(1, max(0, (pmax - math.log(M)) / (6 * math.log(10))))
+    dn, med = disp(ts)
+    if dn is None: return 0.0, 'degenerate_deltas', None
+    reg = min(1, max(0, 1 - 2 * dn))
+    g = min(1.0, reg / gate)
+    per = max(p for p, v in zip(grid, Ps) if v >= .5 * pmax)
+    return 0.5 * sig * g + 0.5 * reg, 'ok', per
+

research/validation/p2real.py

@@ -0,0 +1,52 @@
+# Pillar 2 candidate metrics on real Zeek conn.log: chained session duration per tuple, persistence score,
+# per-host robust z of ln(duration), PCR, failed attempts.
+import sys, math, collections, statistics
+def read(path):
+    fields = None
+    for line in open(path, errors='replace'):
+        line = line.rstrip('\n')
+        if line.startswith('#fields'): fields = line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        yield dict(zip(fields, line.split('\t')))
+num = lambda x: float(x) if x not in ('-', '', '(empty)') else 0.0
+GAP = 3600.0
+def sessions(rows):
+    by4 = collections.defaultdict(list)
+    for r in rows:
+        if r['proto'] != 'tcp': continue
+        by4[(r['id.orig_h'], r['id.orig_p'], r['id.resp_h'], r['id.resp_p'])].append(r)
+    out = []
+    for k, rs in by4.items():
+        rs.sort(key=lambda r: float(r['ts']))
+        cur = None
+        for r in rs:
+            s = float(r['ts']); e = s + num(r['duration'])
+            cont = cur is not None and 'S' not in r['history'] and s - cur['end'] <= GAP
+            if cont:
+                cur['end'] = max(cur['end'], e); cur['ob'] += num(r['orig_bytes']); cur['rb'] += num(r['resp_bytes']); cur['parts'] += 1
+            else:
+                if cur: out.append(cur)
+                cur = dict(key=k, start=s, end=e, ob=num(r['orig_bytes']), rb=num(r['resp_bytes']), parts=1)
+        out.append(cur)
+    return out
+for path in sys.argv[1:]:
+    rows = list(read(path))
+    ss = sessions(rows)
+    fails = collections.Counter((r['id.orig_h'], r['id.resp_h'], r['id.resp_p']) for r in rows if r['proto'] == 'tcp' and r['conn_state'] in ('S0', 'REJ', 'RSTOS0', 'SH'))
+    host_ld = collections.defaultdict(list)
+    for s in ss:
+        d = s['end'] - s['start']
+        if d > 0: host_ld[s['key'][0]].append(math.log(d))
+    best = {}
+    for s in ss:
+        t = (s['key'][0], s['key'][2], s['key'][3]); d = s['end'] - s['start']
+        if t not in best or d > best[t]['end'] - best[t]['start']: best[t] = s
+    ranked = sorted(best.items(), key=lambda kv: -(kv[1]['end'] - kv[1]['start']))
+    print(f"== {path}: tcp sessions={len(ss)} tuples={len(best)}")
+    for t, s in ranked[:5]:
+        d = s['end'] - s['start']; xs = host_ld[t[0]]
+        xs = xs or [0.0]; med = statistics.median(xs); mad = statistics.median(abs(x - med) for x in xs) or 1e-9
+        z = (math.log(d) - med) / (1.4826 * mad) if d > 0 else 0
+        pers = min(1, max(0, math.log(d / 3600) / math.log(24))) if d > 0 else 0
+        tot = s['ob'] + s['rb']; pcr = (s['ob'] - s['rb']) / tot if tot else 0
+        print(f"  {t} L={d:.0f}s ({d/3600:.2f}h) parts={s['parts']} persistence={pers:.3f} host_robust_z={z:.1f} ob={s['ob']:.0f} rb={s['rb']:.0f} PCR={pcr:+.2f} B/h={(tot/(d/3600)) if d>0 else 0:.0f} failed_attempts={fails[t]}")

research/validation/p3hourly.py

@@ -0,0 +1,27 @@
+# Hourly (UTC-aligned) information volume of unique subdomains per (host, base) on real dns.log files.
+import sys, os, collections, math, importlib.util
+here = os.path.dirname(os.path.abspath(__file__))
+sp = importlib.util.spec_from_file_location('p3lib', os.path.join(here, 'p3lib.py')); L = importlib.util.module_from_spec(sp); sp.loader.exec_module(L)
+def read(path):
+    fields = None
+    for line in open(path, errors='replace'):
+        line = line.rstrip('\n')
+        if line.startswith('#fields'): fields = line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        yield dict(zip(fields, line.split('\t')))
+allmax = []
+for path in sys.argv[1:]:
+    b = collections.defaultdict(set)
+    for r in read(path):
+        q = r.get('query')
+        if q in (None, '-', ''): continue
+        base, sub = L.split(q)
+        if sub: b[(r['id.orig_h'], base, int(float(r['ts']) // 3600))].add(sub)
+    best = []
+    for (h, base, hr), subs in b.items():
+        cnt = collections.Counter(c for s in subs for c in L.chars(s))
+        bits = L.H(cnt) * sum(len(L.chars(s)) for s in subs)
+        best.append((bits, h, base, hr, len(subs)))
+    best.sort(reverse=True)
+    print(f"== {path}: (host,base,hour) groups={len(best)}; top3 hourly info bits:")
+    for x in best[:3]: print(f"   bits={x[0]:.0f} host={x[1]} base={x[2]} U={x[4]}")

research/validation/p3lib.py

@@ -0,0 +1,42 @@
+import math, collections
+SLD2 = {'co', 'com', 'net', 'org', 'gov', 'edu', 'ac'}
+def split(q):
+    q = q.rstrip('.').lower()
+    labels = q.split('.')
+    k = 3 if len(labels) >= 3 and len(labels[-1]) == 2 and labels[-2] in SLD2 else 2
+    if len(labels) <= k: return q, ''
+    return '.'.join(labels[-k:]), '.'.join(labels[:-k])
+def chars(s): return [c for c in s if c != '.']
+def H(counter):
+    n = sum(counter.values())
+    return -sum(v / n * math.log2(v / n) for v in counter.values()) if n else 0.0
+def norm(counter):
+    n = sum(counter.values()); return {k: v / n for k, v in counter.items()} if n else {}
+def jsd(p, q):
+    keys = sorted(set(p) | set(q)); s = 0.0
+    for k in keys:
+        a = p.get(k, 0.0); b = q.get(k, 0.0); m = (a + b) / 2
+        if a > 0: s += 0.5 * a * math.log2(a / m)
+        if b > 0: s += 0.5 * b * math.log2(b / m)
+    return s
+def groups_from(pairs):
+    g = collections.defaultdict(set)
+    for host, q in pairs:
+        base, sub = split(q)
+        if sub: g[(host, base)].add(sub)
+    return g
+def features(g):
+    dists = {k: norm(collections.Counter(c for s in subs for c in chars(s))) for k, subs in g.items()}
+    keys = sorted(dists); N = len(keys)
+    tot = collections.defaultdict(float)
+    for k in keys:
+        for c, v in dists[k].items(): tot[c] += v
+    out = {}
+    for k in keys:
+        subs = g[k]
+        ref = {c: (tot[c] - dists[k].get(c, 0.0)) / (N - 1) for c in tot} if N > 1 else {}
+        ref = {c: v for c, v in ref.items() if v > 1e-15}
+        cnt = collections.Counter(c for s in subs for c in chars(s))
+        out[k] = dict(U=len(subs), H=H(cnt), JSD=jsd(dists[k], ref) if ref else 0.0,
+                      L=sum(len(s) for s in subs) / len(subs), maxlab=max(max(len(x) for x in s.split('.')) for s in subs))
+    return out

research/validation/p3real.py

@@ -0,0 +1,18 @@
+import sys, os, importlib.util
+here = os.path.dirname(os.path.abspath(__file__))
+sp = importlib.util.spec_from_file_location('p3lib', os.path.join(here, 'p3lib.py')); L = importlib.util.module_from_spec(sp); sp.loader.exec_module(L)
+def read(path):
+    fields = None
+    for line in open(path, errors='replace'):
+        line = line.rstrip('\n')
+        if line.startswith('#fields'): fields = line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        yield dict(zip(fields, line.split('\t')))
+for path in sys.argv[1:]:
+    pairs = [(r['id.orig_h'], r['query']) for r in read(path) if r.get('query') not in (None, '-', '')]
+    g = L.groups_from(pairs); f = L.features(g)
+    top = sorted(f.items(), key=lambda kv: -kv[1]['U'])[:6]
+    print(f"== {path}: queries={len(pairs)} groups={len(f)}")
+    for k, v in top: print(f"  {k} U={v['U']} H={v['H']:.2f} JSD={v['JSD']:.3f} meanlen={v['L']:.1f} maxlabel={v['maxlab']}")
+    hj = sorted(f.items(), key=lambda kv: -kv[1]['JSD'])[:3]
+    for k, v in hj: print(f"  topJSD {k} U={v['U']} H={v['H']:.2f} JSD={v['JSD']:.3f} meanlen={v['L']:.1f}")

research/validation/p3sim.py

@@ -0,0 +1,37 @@
+# Inject real tunnel qnames (Mendeley DTQBC, CC BY 4.0) at various volumes into real benign DNS (CTU-Normal-20 dns.log
+# regenerated with Zeek 9.0.1); compute per-(host, base) features. 300 seeds per (class, U).
+import sys, os, csv, random, statistics, math, importlib.util
+here = os.path.dirname(os.path.abspath(__file__))
+sp = importlib.util.spec_from_file_location('p3lib', os.path.join(here, 'p3lib.py')); L = importlib.util.module_from_spec(sp); sp.loader.exec_module(L)
+def read(path):
+    fields = None
+    for line in open(path, errors='replace'):
+        line = line.rstrip('\n')
+        if line.startswith('#fields'): fields = line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        yield dict(zip(fields, line.split('\t')))
+benign_logs = sys.argv[1].split(',')
+mend = sys.argv[2]; seeds = int(sys.argv[3])
+rows = list(csv.reader(open(mend, encoding='utf-8', errors='replace')))[1:]
+cls = {c: [r[1] for r in rows if r[0] == c] for c in '1234'}
+names = {'1': 'dns2tcp', '2': 'dnscapy', '3': 'iodine', '4': 'tuns'}
+def info_bits(f): return f['H'] * f['L'] * f['U']
+for bl in benign_logs:
+    base_pairs = [(r['id.orig_h'], r['query']) for r in read(bl) if r.get('query') not in (None, '-', '')]
+    fb = L.features(L.groups_from(base_pairs))
+    ib = sorted((info_bits(v), k, v) for k, v in fb.items())
+    print(f"== benign {bl}: groups={len(fb)}; top info_bits:")
+    for i, k, v in ib[-4:]: print(f"   {k} U={v['U']} H={v['H']:.2f} meanlen={v['L']:.1f} JSD={v['JSD']:.3f} info_bits={i:.0f}")
+    for c in '124':
+        for U in (30, 60, 100, 300):
+            I, J, Hs, Ls = [], [], [], []
+            for s in range(1, seeds + 1):
+                r = random.Random(s * 1000 + U)
+                pool = list(dict.fromkeys(cls[c])); sample = r.sample(pool, min(U, len(pool)))
+                pairs = base_pairs + [('10.66.66.66', q) for q in sample]
+                f = L.features(L.groups_from(pairs))
+                tun = [v for k, v in f.items() if k[0] == '10.66.66.66']
+                v = max(tun, key=lambda x: x['U'])
+                I.append(info_bits(v)); J.append(v['JSD']); Hs.append(v['H']); Ls.append(v['L'])
+            q = lambda xs: f"min={min(xs):.3g} med={statistics.median(xs):.3g} max={max(xs):.3g}"
+            print(f"  {names[c]:8s} U={U:4d} info_bits {q(I)} | JSD {q(J)} | H {q(Hs)} | meanlen {q(Ls)}")

research/validation/zeekstat.py

@@ -0,0 +1,38 @@
+import sys, collections
+def read(path):
+    fields=None; rows=[]
+    for line in open(path, errors='replace'):
+        line=line.rstrip('\n')
+        if line.startswith('#fields'): fields=line.split('\t')[1:]; continue
+        if line.startswith('#') or not line: continue
+        p=line.split('\t')
+        rows.append(dict(zip(fields,p)))
+    return rows
+def f(x):
+    try: return float(x)
+    except: return None
+for path in sys.argv[1:]:
+    rows=read(path)
+    tup=collections.defaultdict(list)
+    for r in rows:
+        if r.get('proto')!='tcp': continue
+        tup[(r['id.orig_h'],r['id.resp_h'],r['id.resp_p'])].append(r)
+    out=[]
+    for k,rs in tup.items():
+        ivs=sorted((f(r['ts']), f(r['ts'])+(f(r['duration']) or 0)) for r in rs)
+        # union length
+        tot=0; cs,ce=ivs[0]
+        for s,e in ivs[1:]:
+            if s<=ce: ce=max(ce,e)
+            else: tot+=ce-cs; cs,ce=s,e
+        tot+=ce-cs
+        mx=max((f(r['duration']) or 0) for r in rs)
+        ob=sum(int(r['orig_bytes']) if r['orig_bytes'].isdigit() else 0 for r in rs)
+        rb=sum(int(r['resp_bytes']) if r['resp_bytes'].isdigit() else 0 for r in rs)
+        states=collections.Counter(r['conn_state'] for r in rs)
+        out.append((tot,mx,len(rs),k,ob,rb,dict(states.most_common(3))))
+    out.sort(reverse=True)
+    allts=[f(r['ts']) for r in rows]
+    print('==',path, 'tcp tuples',len(out),'span_h %.2f'%((max(allts)-min(allts))/3600))
+    for o in out[:8]:
+        print('  union=%.0fs max=%.0fs n=%d %s ob=%d rb=%d %s'%o)