#!/usr/bin/env python3
"""Turns the three instruments into the numbers the article may use.

A file is scored only when the fetch returned 200 or 206 and the parser recorded no error, which is
the denominator rule the pre-registration fixed; the excluded count travels beside every share.

A share over indicator data has an exact bootstrap: resampling n indicators with replacement is the
same draw as Binomial(n, p)/n. A difference between two independent groups is the difference of two
such draws, and a paired difference is a multinomial draw over the four cells. Every interval here
is 10 000 real resamples.
"""
import json, os
from collections import defaultdict, Counter
import numpy as np

RUN = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
D = os.path.join(RUN, 'data')
BOOT = 10000
RNG = np.random.default_rng(20260910)

# Fixed here, before the full frame was read, and named in the article's method: a make is a phone
# make when the devices it sells under that name are predominantly smartphones. The list is the
# study's, not the data's, and the per-make table beside every family figure lets a reader redo the
# grouping.
PHONE_MAKES = {'Apple', 'samsung', 'SAMSUNG', 'Xiaomi', 'Google', 'HUAWEI', 'motorola', 'Motorola',
               'Nokia', 'LG Electronics', 'LGE', 'OnePlus', 'OPPO', 'realme', 'ZTE', 'vivo',
               'Sony Ericsson', 'HTC', 'BlackBerry', 'Wiko', 'alcatel', 'ASUS', 'Lenovo'}


def load(name):
    p = os.path.join(D, name)
    if not os.path.exists(p):
        return []
    with open(p) as f:
        return [json.loads(l) for l in f if l.strip()]


def boot_share(flags, boot=BOOT):
    x = np.asarray(flags, dtype=float)
    if len(x) == 0:
        return None
    draws = RNG.binomial(len(x), x.mean(), size=boot) / len(x)
    return [round(float(np.percentile(draws, 2.5)), 4), round(float(np.percentile(draws, 97.5)), 4)]


def boot_two_groups(a, b, boot=BOOT):
    """95 per cent interval on the difference between two independent groups."""
    a = np.asarray(a, float); b = np.asarray(b, float)
    if len(a) == 0 or len(b) == 0:
        return None
    da = RNG.binomial(len(a), a.mean(), size=boot) / len(a)
    db = RNG.binomial(len(b), b.mean(), size=boot) / len(b)
    d = da - db
    return [round(float(np.percentile(d, 2.5)), 4), round(float(np.percentile(d, 97.5)), 4)]


def boot_paired(table, boot=BOOT):
    """95 per cent interval on a - b for two indicators on the same files.

    table = {'both': x, 'a_only': x, 'b_only': x, 'neither': x}. Resampling the files is a
    multinomial draw with those four probabilities.
    """
    counts = np.array([table['both'], table['a_only'], table['b_only'], table['neither']], float)
    n = counts.sum()
    if n == 0:
        return None
    draws = RNG.multinomial(int(n), counts / n, size=boot)
    d = (draws[:, 0] + draws[:, 1] - draws[:, 0] - draws[:, 2]) / n
    return [round(float(np.percentile(d, 2.5)), 4), round(float(np.percentile(d, 97.5)), 4)]


def share(num, den):
    return round(num / den, 4) if den else None


def sanitize(key):
    out = ''.join(c if (c.isalnum() or c in '_-') else '_' for c in key)
    return out if not out[:1].isdigit() else 'm_' + out


def scored(rows):
    return [r for r in rows if r.get('http_status') in (200, 206) and not r.get('parse_error')]


def file_serial_any(r):
    """What a reader who opens the file itself finds, without the maker's own note."""
    return bool(r.get('body_serial_tag_present') or r.get('xmp_serial_present'))


def has_gps(r):
    return bool(r.get('gps_coord_tag_present') and not r.get('gps_coord_all_zero'))


# the whole-file reading by title, for comparing it with the wider exiftool pass
whole_by_title = {r['title']: (r.get('whole') or {})
                  for r in [json.loads(l) for l in open(os.path.join(D, 'validation.jsonl')) if l.strip()]}


units = load('units.jsonl')
validation = load('validation.jsonl')
frame_summary = json.load(open(os.path.join(D, 'frame-summary.json')))
cross = {r['title']: r for r in json.load(open(os.path.join(D, 'crosscheck.json')))['rows']}
sc = scored(units)

out = {
    'run_date': '2026-09-10',
    'head_window_bytes': 131072,
    'generated_at': '2026-09-10',
    'source': 'Wikimedia Commons file namespace, random draw, JPEG files only',
    'boot_resamples': BOOT,
}

F = {
    'drawn': frame_summary.get('drawn'),
    'accepted': frame_summary.get('accepted'),
    'acceptance_rate': frame_summary.get('acceptance_rate'),
    'rejections': frame_summary.get('rejections'),
    'rows': len(units), 'scored': len(sc), 'excluded': len(units) - len(sc),
}
for name, key, fn in (
    ('make', 'make_present', lambda r: bool(r.get('make_present'))),
    ('model', 'model_present', lambda r: bool(r.get('model_present'))),
    ('exif', 'has_exif', lambda r: bool(r.get('has_exif'))),
    ('gps_tag', 'gps_coord_tag_present', lambda r: bool(r.get('gps_coord_tag_present'))),
    ('gps_coord', 'gps_coord_coord', lambda r: has_gps(r)),
    ('gps_tag_all_zero', 'gps_tag_all_zero', lambda r: bool(r.get('gps_coord_tag_present')) and bool(r.get('gps_coord_all_zero'))),
    ('datetime', 'datetime_original_present', lambda r: bool(r.get('datetime_original_present'))),
    ('software', 'software_present', lambda r: bool(r.get('software_present'))),
    ('serial_body', 'body_serial_tag_present', lambda r: bool(r.get('body_serial_tag_present'))),
    ('serial_lens', 'lens_serial_tag_present', lambda r: bool(r.get('lens_serial_tag_present'))),
    ('serial_xmp', 'xmp_serial_present', lambda r: bool(r.get('xmp_serial_present'))),
    ('serial_any', 'serial_any', lambda r: file_serial_any(r)),
):
    flags = [1 if fn(r) else 0 for r in sc]
    F['with_' + name] = sum(flags)
    F['share_' + name] = share(sum(flags), len(sc))
    if name in ('make', 'gps_coord', 'serial_body', 'serial_xmp', 'serial_any', 'exif'):
        F['ci_' + name] = boot_share(flags)
out['frame'] = F

# --- the platform's own reading, over the same frame ------------------------------
paired = [r for r in sc if r['title'] in cross]
table = {'both': 0, 'a_only': 0, 'b_only': 0, 'neither': 0}
for r in paired:
    a = bool(cross[r['title']].get('serial'))
    b = file_serial_any(r)
    table['both' if (a and b) else 'a_only' if a else 'b_only' if b else 'neither'] += 1
n = len(paired)
out['paired'] = {
    'n': n, 'a_name': 'platform', 'b_name': 'file',
    'a_with': table['both'] + table['a_only'], 'b_with': table['both'] + table['b_only'],
    'a_share': share(table['both'] + table['a_only'], n),
    'b_share': share(table['both'] + table['b_only'], n),
    'difference_points': share(table['a_only'] - table['b_only'], n),
    'file_minus_platform_points': share(table['b_only'] - table['a_only'], n),
    # the interval above is for platform minus file; this one is for file minus platform
    'ci_file_minus_platform': [-c for c in reversed(boot_paired(table) or [])],
    'ci_difference': boot_paired(table),
    **table,
}
out['platform'] = {
    'scored': n,
    'with_serial': table['both'] + table['a_only'],
    'share_serial': share(table['both'] + table['a_only'], n),
    'with_make': sum(1 for r in paired if cross[r['title']].get('make')),
    'share_make': share(sum(1 for r in paired if cross[r['title']].get('make')), n),
    'with_gps': sum(1 for r in paired if cross[r['title']].get('gps')),
    'share_gps': share(sum(1 for r in paired if cross[r['title']].get('gps')), n),
    'with_datetime': sum(1 for r in paired if cross[r['title']].get('datetime_original')),
    'share_datetime': share(sum(1 for r in paired if cross[r['title']].get('datetime_original')), n),
    'with_software': sum(1 for r in paired if cross[r['title']].get('software')),
    'share_software': share(sum(1 for r in paired if cross[r['title']].get('software')), n),
}

# --- the whole-file subsample, with exiftool as referee ---------------------------
val = [r for r in validation if (r.get('exiftool') or {}).get('groups') is not None]


def ex_groups(r):
    return set(g.split(':')[0] for g in (r.get('exiftool') or {}).get('serial_where', []))


whole = lambda r: (r.get('whole') or {})
V = {
    'n': len(val),
    'capped_by_size': sum(1 for r in validation if r.get('capped')),
    'our_body': sum(1 for r in val if whole(r).get('body_serial_tag_present')),
    'our_xmp': sum(1 for r in val if whole(r).get('xmp_serial_present')),
    'our_any': sum(1 for r in val if whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present')),
    'exiftool_any': sum(1 for r in val if ex_groups(r)),
    'exiftool_exif': sum(1 for r in val if 'EXIF' in ex_groups(r)),
    'exiftool_xmp': sum(1 for r in val if 'XMP' in ex_groups(r)),
    'exiftool_makernote': sum(1 for r in val if 'MakerNotes' in ex_groups(r)),
    'in_makernote_only': sum(1 for r in val if ex_groups(r) == {'MakerNotes'}),
    'in_xmp_only': sum(1 for r in val if ex_groups(r) == {'XMP'}),
    'in_exif_only': sum(1 for r in val if ex_groups(r) == {'EXIF'}),
    'missed_by_our_reading': sum(1 for r in val if ex_groups(r) and not (whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present'))),
    'our_reading_without_any': sum(1 for r in val if not ex_groups(r) and (whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present'))),
}
for k in ('our_any', 'exiftool_any', 'exiftool_makernote', 'missed_by_our_reading'):
    V['share_' + k] = share(V[k], len(val))
    V['ci_' + k] = boot_share([V[k + '_flags'][i] for i in range(len(val))]) if False else None
V['ci_exiftool_any'] = boot_share([1 if ex_groups(r) else 0 for r in val])
V['ci_our_any'] = boot_share([1 if (whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present')) else 0 for r in val])
V['ci_makernote'] = boot_share([1 if 'MakerNotes' in ex_groups(r) else 0 for r in val])
pt = {'both': 0, 'a_only': 0, 'b_only': 0, 'neither': 0}
for r in val:
    a = bool(ex_groups(r))
    b = bool(whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present'))
    pt['both' if (a and b) else 'a_only' if a else 'b_only' if b else 'neither'] += 1
V['difference_points'] = share(pt['a_only'] - pt['b_only'], len(val))
V['ci_difference'] = boot_paired(pt)

# the head window, measured against the whole file on the same files
head_by_title = {r['title']: r for r in load('units.jsonl')}
head_missed_xmp = sum(1 for r in val if whole(r).get('xmp_serial_present')
                      and not (head_by_title.get(r['title']) or {}).get('xmp_serial_present'))
V['head_window_missed_xmp'] = head_missed_xmp
# The claim is about serial numbers the whole file holds, so it is the union of the two fields this
# study reads, not the XMP field alone: a body serial written past the window counts too.
def _our_serial(r):
    return bool(whole(r).get('body_serial_tag_present') or whole(r).get('xmp_serial_present'))
def _head_serial(r):
    h = head_by_title.get(r['title']) or {}
    return bool(h.get('body_serial_tag_present') or h.get('xmp_serial_present'))
V['head_window_missed'] = sum(1 for r in val if _our_serial(r) and not _head_serial(r))
V['head_window_missed_share'] = share(V['head_window_missed'], len(val))
# A bootstrap of an all-zero sample returns zero at both ends, which states no boundary at all.
# An every-file check needs the exact one-sided bound instead: how large a miss rate a clean run of
# this many files still permits.
_k = V['head_window_missed']
try:
    from scipy.stats import beta as _beta
    _hi = float(_beta.ppf(0.975, _k + 1, len(val) - _k)) if _k < len(val) else 1.0
except Exception:
    _hi = 1 - 0.05 ** (1.0 / len(val))
V['ci_head_window_missed'] = [share(_k, len(val)), _hi]
V['ci_head_window_missed_method'] = (
    'an exact one-sided bound: a clean run of %d files still permits a miss rate this high' % len(val))
V['head_window_missed_make'] = sum(1 for r in val if (whole(r).get('make') or '')
                                   and not (head_by_title.get(r['title']) or {}).get('make'))
V['whole_with_xmp'] = sum(1 for r in val if whole(r).get('xmp_serial_present'))
# which group the serial sits in, per maker, on the whole-file subsample
per_make_groups = defaultdict(lambda: defaultdict(int))
for r in val:
    m = (whole(r).get('make') or '').strip() or '(no make)'
    m = sanitize(m)
    g = ex_groups(r)
    per_make_groups[m]['n'] += 1
    for grp in ('EXIF', 'XMP', 'MakerNotes'):
        if grp in g:
            per_make_groups[m][grp.lower()] += 1
    if g:
        per_make_groups[m]['any'] += 1
V['by_make'] = {k: dict(v) for k, v in sorted(per_make_groups.items(), key=lambda kv: -kv[1]['n'])}
out['validation'] = V

# --- the wider exiftool reading, run after the second method exposed the gap ----------------
# The first pass asked exiftool for SerialNumber and BodySerialNumber only. A second, independently
# derived method asked for InternalSerialNumber and CameraSerialNumber as well and found a serial in
# far more files. This block is that reading, on the same files and by the same program, and it is
# what the article binds to: the narrower pass is kept in the ledger as the row it supersedes.
wide = [r for r in load('validation-wide.jsonl') if not r.get('error')]
W = {'n': len(wide)}
if wide:
    def where(r):
        # A lens serial (0xA435) identifies the lens, not the camera. It is the trap this run
        # already recorded as D-06, and it reappears here because this pass asks exiftool for
        # more tags than the narrow one did, the lens among them.
        return [s for s in (r.get('serial_where') or [])
                if s and not s.split(':')[-1] == 'LensSerialNumber']
    def any_serial(r): return bool(where(r))
    def named(r): return [s for s in where(r) if 'internal' not in s.lower()]
    def internal(r): return [s for s in where(r) if 'internal' in s.lower()]
    def in_group(r, name):
        # the tag is named with its group, so the group of a serial is read from the tag itself
        return any(s.split(':')[0] == name for s in where(r))
    W['any'] = sum(1 for r in wide if any_serial(r))
    W['share_any'] = share(W['any'], len(wide))
    W['ci_any'] = boot_share([1 if any_serial(r) else 0 for r in wide])
    W['named_only'] = sum(1 for r in wide if named(r) and not internal(r))
    W['share_named'] = share(sum(1 for r in wide if named(r)), len(wide))
    W['internal_only'] = sum(1 for r in wide if internal(r) and not named(r))
    W['share_internal'] = share(W['internal_only'], len(wide))
    W['both'] = sum(1 for r in wide if named(r) and internal(r))
    W['with_make'] = sum(1 for r in wide if r.get('has_make'))
    W['in_makernote'] = sum(1 for r in wide if in_group(r, 'MakerNotes'))
    W['in_exif'] = sum(1 for r in wide if in_group(r, 'EXIF'))
    W['in_xmp'] = sum(1 for r in wide if in_group(r, 'XMP'))
    W['makernote_only'] = sum(1 for r in wide if in_group(r, 'MakerNotes')
                              and not in_group(r, 'EXIF') and not in_group(r, 'XMP'))
    W['exif_only'] = sum(1 for r in wide if in_group(r, 'EXIF')
                         and not in_group(r, 'MakerNotes') and not in_group(r, 'XMP'))
    W['xmp_only'] = sum(1 for r in wide if in_group(r, 'XMP')
                        and not in_group(r, 'EXIF') and not in_group(r, 'MakerNotes'))
    W['missed_by_our_reading'] = sum(
        1 for r in wide if any_serial(r) and not (whole_by_title.get(r['title'], {}).get('body_serial_tag_present')
                                                  or whole_by_title.get(r['title'], {}).get('xmp_serial_present')))
    # the same three-way comparison against this study's own reading, on the same files
    pt2 = {'both': 0, 'a_only': 0, 'b_only': 0, 'neither': 0}
    for r in wide:
        a = any_serial(r)
        w = whole_by_title.get(r['title'], {})
        b = bool(w.get('body_serial_tag_present') or w.get('xmp_serial_present'))
        pt2['both' if (a and b) else 'a_only' if a else 'b_only' if b else 'neither'] += 1
    W['our_reading_n'] = pt2['both'] + pt2['b_only']
    W['share_our_reading'] = share(W['our_reading_n'], len(wide))
    W['difference_points'] = share(pt2['a_only'] - pt2['b_only'], len(wide))
    W['ci_difference'] = boot_paired(pt2)
    W['not_seen_by_either'] = pt2['neither']
    W.update({f'cells_{k}': v for k, v in pt2.items()})
out['validation_wide'] = W

# --- the coordinate on the page rather than in the file --------------------------------------
pc = {f['title']: f.get('page_coordinate') for f in
      [json.loads(l) for l in open(os.path.join(D, 'frame.jsonl')) if l.strip()]}
flags = [1 if pc.get(r['title']) else 0 for r in sc]
# A page coordinate and a file coordinate are different things on the same file, so the gap between
# them is a paired difference and gets the paired interval.
_page_cells = {'both': 0, 'a_only': 0, 'b_only': 0, 'neither': 0}
for r in sc:
    a = bool(pc.get(r['title']))
    b = has_gps(r)
    _page_cells['both' if (a and b) else 'a_only' if a else 'b_only' if b else 'neither'] += 1
out['page_coordinate'] = {
    'n': len(sc), 'with_page_coordinate': sum(flags), 'share': share(sum(flags), len(sc)),
    'ci': boot_share(flags),
    'with_exif_coordinate': F['with_gps_coord'], 'share_exif': F['share_gps_coord'],
    'difference_points': share(_page_cells['a_only'] - _page_cells['b_only'], len(sc)),
    'ci_difference': boot_paired(_page_cells),
    'cells': _page_cells,
}

# --- per make --------------------------------------------------------------------
by_make = defaultdict(lambda: Counter())
for r in sc:
    m = (r.get('make') or '').strip()
    if not m:
        continue
    c = by_make[m]
    c['n'] += 1
    c['body'] += bool(r.get('body_serial_tag_present'))
    c['xmp'] += bool(r.get('xmp_serial_present'))
    c['lens'] += bool(r.get('lens_serial_tag_present'))
    c['any'] += file_serial_any(r)
    c['gps'] += has_gps(r)
    c['model'] += bool(r.get('model_present'))

makes = {}
for m, c in by_make.items():
    makes[sanitize(m)] = {
        'label': m, 'n': c['n'], 'files': c['n'],
        'share_of_frame': share(c['n'], len(sc)),
        'with_serial_body': c['body'], 'share_serial_body': share(c['body'], c['n']),
        'with_serial_xmp': c['xmp'], 'share_serial_xmp': share(c['xmp'], c['n']),
        'with_serial_lens': c['lens'], 'share_serial_lens': share(c['lens'], c['n']),
        'with_serial_any': c['any'], 'share_serial_any': share(c['any'], c['n']),
        'with_gps': c['gps'], 'share_gps': share(c['gps'], c['n']),
        'with_model': c['model'],
        'is_phone': m in PHONE_MAKES,
    }
# A zero cell still needs a boundary. The per-make table shows one: the bare NIKON name carries a
# serial in none of its files, and the reader cannot see from a zero how wide that statement is.
_nk = makes.get('NIKON')
if _nk:
    _n, _k = _nk['n'], _nk['with_serial_any']
    try:
        from scipy.stats import beta as _b2
        _hi = float(_b2.ppf(0.975, _k + 1, _n - _k)) if _k < _n else 1.0
    except Exception:
        _hi = 1 - 0.05 ** (1.0 / _n)
    out['zero_cell'] = {
        'make': 'NIKON', 'n': _n, 'with_serial': _k,
        'share': share(_k, _n), 'upper': _hi,
    }

out['by_make'] = makes
out['makes_by_size'] = [k for k, v in sorted(makes.items(), key=lambda kv: -kv[1]['n'])]

# --- camera against phone ---------------------------------------------------------
groups = {'camera': [], 'phone': []}
for r in sc:
    m = (r.get('make') or '').strip()
    if not m:
        continue
    groups['phone' if m in PHONE_MAKES else 'camera'].append(r)

fam = {}
for name, rows in groups.items():
    s = [1 if file_serial_any(r) else 0 for r in rows]
    g = [1 if has_gps(r) else 0 for r in rows]
    ex = [1 if r.get('has_exif') else 0 for r in rows]
    mk = [1 if r.get('make_present') else 0 for r in rows]
    dt = [1 if r.get('datetime_original_present') else 0 for r in rows]
    sw = [1 if r.get('software_present') else 0 for r in rows]
    fam[name] = {
        'n': len(rows),
        'with_datetime': sum(dt), 'share_datetime': share(sum(dt), len(rows)),
        'with_software': sum(sw), 'share_software': share(sum(sw), len(rows)),
        'with_exif': sum(ex), 'share_exif': share(sum(ex), len(rows)),
        'with_make': sum(mk), 'share_make': share(sum(mk), len(rows)),
        'with_serial': sum(s), 'share_serial': share(sum(s), len(rows)),
        'ci_serial': boot_share(s),
        'with_gps': sum(g), 'share_gps': share(sum(g), len(rows)),
        'ci_gps': boot_share(g),
    }
fam['serial_difference_points'] = round(fam['camera']['share_serial'] - fam['phone']['share_serial'], 4)
fam['ci_serial_difference'] = boot_two_groups(
    [1 if file_serial_any(r) else 0 for r in groups['camera']],
    [1 if file_serial_any(r) else 0 for r in groups['phone']])
fam['gps_difference_points'] = round(fam['phone']['share_gps'] - fam['camera']['share_gps'], 4)
fam['ci_gps_difference'] = boot_two_groups(
    [1 if has_gps(r) else 0 for r in groups['phone']],
    [1 if has_gps(r) else 0 for r in groups['camera']])
out['families'] = fam

# --- per model, above a stated size ------------------------------------------------
# A table of 745 models is a dataset, not a table, so the article carries the models with at least
# twenty files, and the threshold is stated beside it rather than chosen to fit the answer.
MODEL_MIN = 20
by_model = defaultdict(lambda: Counter())
for r in sc:
    m = (r.get('model') or '').strip()
    if not m:
        continue
    c = by_model[m]
    c['n'] += 1
    c['any'] += file_serial_any(r)
    c['body'] += bool(r.get('body_serial_tag_present'))
    c['gps'] += has_gps(r)
models = {}
for m, c in by_model.items():
    models[sanitize(m)] = {
        'label': m, 'n': c['n'],
        'with_serial': c['any'], 'share_serial': share(c['any'], c['n']),
        'with_serial_body': c['body'],
        'with_gps': c['gps'], 'share_gps': share(c['gps'], c['n']),
        'make': sanitize((next((r.get('make') or '' for r in sc if (r.get('model') or '').strip() == m), '') or '').strip()),
    }
out['model_min'] = MODEL_MIN
out['models_at_or_above_min'] = sorted([k for k, v in models.items() if v['n'] >= MODEL_MIN],
                                       key=lambda k: -models[k]['n'])
out['models_distinct'] = len(models)
out['by_model'] = models

# --- concentration ------------------------------------------------------------------
carriers = sorted([(m, c['any']) for m, c in by_make.items() if c['any'] > 0], key=lambda kv: -kv[1])
total_any = sum(c for _, c in carriers)
top3 = sum(c for _, c in carriers[:3])
big = {m: c for m, c in by_make.items() if c['n'] >= 20}
out['concentration'] = {
    'distinct_makes_with_a_serial': len(carriers),
    'distinct_makes_scored': len(by_make),
    'files_with_serial': total_any,
    'top3_makes': [(m, c) for m, c in carriers[:3]],
    'top3_files': top3,
    'top3_share': share(top3, total_any),
    'makes_with_20_or_more': len(big),
    'makes_at_or_below_one_percent': sorted([m for m, c in big.items() if c['any'] / c['n'] <= 0.01]),
    'makes_at_or_below_one_percent_count': len([m for m, c in big.items() if c['any'] / c['n'] <= 0.01]),
    'overall_serial_share': F['share_serial_any'],
    'max_gap_points': max((round(abs(c['any'] / c['n'] - F['share_serial_any']), 4) for c in big.values()), default=None),
    'largest_make_without_a_serial': sorted([(m, c['n']) for m, c in by_make.items() if c['any'] == 0], key=lambda kv: -kv[1])[:6],
}

with open(os.path.join(D, 'analysis.json'), 'w') as f:
    json.dump(out, f, indent=2)
print(json.dumps({'frame': {k: v for k, v in F.items() if k.startswith(('share', 'scored', 'with_serial', 'ci_'))},
                  'paired': out['paired'], 'families': fam, 'concentration': out['concentration'],
                  'validation': {k: v for k, v in V.items() if k != 'by_make'}},
                 indent=2))