"""Between-group separation B and within-group early/late discrepancy E (two candidates, two weightings).
This ratio B/E is a temporal separation diagnostic. It is NOT HiFi's fidelity F = L/R, and the
names here never say "fidelity"; :mod:`.archived_names` maps them to the archived
temporal-fidelity-Pareto (PF) column names for regression checks only.
Definitions (per paired early/late dataset; frames are the unit)
----------------------------------------------------------------
* Joint set ``G_joint = {g : n_g^E >= 1 and n_g^L >= 1}`` (complete visits). Every candidate and
weighting uses this set; the same weights serve both halves, numerator and denominator.
* ``w_visit = (n_E + n_L) / T`` with ``T`` the joint total, ``w_equal = 1/|G_joint|``.
* ``B_tau = 1/2 sum_{g,h} w_g w_h (|mu_g^E - mu_h^E| + |mu_g^L - mu_h^L|)`` over ordered pairs with
the zero diagonal; ``E_tau = sum_g w_g |mu_g^E - mu_g^L|``.
* ``B_P``/``E_P``: the same with the exact empirical W1 distance of complete-visit durations.
Production arithmetic (a port of PF's canonical float algorithm, same operation order)
---------------------------------------------------------------------------------------
``S_w(v) = 2 sum_b (m_b/T) A_b`` with ``A_1 = 0``, ``A_b = A_{b-1} + (M_{<b}/T)(v_b - v_{b-1})`` over
exact-tie blocks of the stably sorted values; every sum is a sequential ascending sum (``np.cumsum``,
never pairwise ``np.sum``). ``B_tau`` uses float means ``s/n``; ``B_P`` integrates ``S_w`` of float
CDF ordinates ``C/n`` over the union support in increasing duration. ``E`` is the sequential sum of
the ascending-sorted per-group terms ``w_g * Delta_g`` with ``Delta_g = N_g / (n_E n_L)`` and EXACT
integer numerators ``N_tau = |s_E n_L - s_L n_E|`` and ``N_W1 = sum_i |C_E n_L - C_L n_E| (u_{i+1}-u_i)``.
A guard from observed maxima runs BEFORE any product; if a bound could exceed int64 the same
algorithm runs on Python integers (``path='python_int'``) with exact, unclamped numerators. The
guard bound ``GUARD_LIMIT = 2**62`` may only be lowered (``guard_limit``), never raised.
Statuses come from the exact integer identities (empty joint set, B = 0, E = 0), never from a
float magnitude threshold; the one float-derived status, ``numerical_indeterminate``, marks
production floats that are non-finite or contradict those identities. Such points go to the
non-finite register; the Level-1 enclosure does not repair them.
Arithmetic certification (labelled ARITHMETIC, not statistical confidence)
---------------------------------------------------------------------------
Level 1 (:func:`level1`): binary64 interval enclosures of exact B and E from exact rational
blocks with one-ulp outward rounding (port of PF ``exact_verify.level1``); orderings and ties
come from exact integer cross-multiplication, so a float-sort disagreement is repaired instead of
aborting. Level 2b (:func:`exact_components`): exact ``Fraction`` values under a budget.
"""
from __future__ import annotations
from fractions import Fraction
from functools import cmp_to_key
import math
import time
import numpy as np
from .runs import BRIDGE_FIELDS, half_stats
STATUS = ('finite', 'undefined_no_between_separation_and_no_window_difference',
'zero_denominator_positive_separation', 'unavailable_empty_joint_set', 'numerical_indeterminate')
FINITE, UNDEFINED_00, DIVERGENT, UNAVAILABLE, INDETERMINATE = range(5)
CANDIDATES = ('tau', 'P')
WEIGHTINGS = ('visit', 'equal')
COMBOS = (('tau', 'visit'), ('tau', 'equal'), ('P', 'visit'), ('P', 'equal'))
GUARD_LIMIT = 2 ** 62
FLOAT_FIELDS = ('sep_tau_early_visit', 'sep_tau_late_visit', 'sep_tau_early_equal', 'sep_tau_late_equal',
'E_tau_visit', 'E_tau_equal',
'sep_P_early_visit', 'sep_P_late_visit', 'sep_P_early_equal', 'sep_P_late_equal',
'E_P_visit', 'E_P_equal')
INT_FIELDS = ('G_joint', 'visit_mass_total', 'tau_B_zero_early', 'tau_B_zero_late', 'tau_E_zero',
'P_B_zero_early', 'P_B_zero_late', 'P_E_zero', 'support_union_early', 'support_union_late',
'max_duration_joint_early', 'max_duration_joint_late', 'valid')
RELIABILITY_KIND = 'ARITHMETIC: binary64 rounding enclosure/error of the exact value; not statistical confidence'
_INT64_MAX = 2 ** 63 - 1
_TWO53 = 2 ** 53
[docs]
class BudgetExhausted(RuntimeError):
"""Exact (Level 2b) evaluation stopped by its recorded budget; the relation stays unresolved."""
# ------------------------------------------------------------------------------- exact integers
def _obj(values):
out = np.empty(len(values), dtype=object)
out[:] = list(values)
return out
def _exact_ints(x, name):
"""int64 array, or an object array of Python ints when a value does not fit int64."""
a = np.asarray(x)
if a.ndim != 1:
raise ValueError(f'{name} must be 1-D')
if a.dtype.kind in 'iu':
if a.dtype == np.uint64 and a.size and int(a.max()) > _INT64_MAX:
return _obj(int(v) for v in a.tolist())
return a.astype(np.int64)
if a.dtype == object:
vals = []
for v in a.tolist():
if isinstance(v, (bool, np.bool_)) or not isinstance(v, (int, np.integer)):
raise TypeError(f'{name} must hold exact integers, found {v!r}')
vals.append(int(v))
if all(-_INT64_MAX - 1 <= v <= _INT64_MAX for v in vals):
return np.asarray(vals, np.int64)
return _obj(vals)
raise TypeError(f'{name} must hold exact integers, got dtype {a.dtype}')
def _as_python_ints(a):
return a if a.dtype == object else _obj(a.tolist())
def _to_f64(a):
"""binary64 conversion (round to nearest even), identical for int64 arrays and Python ints."""
a = np.asarray(a)
if a.dtype == object:
return np.asarray([float(v) for v in a.ravel().tolist()], np.float64).reshape(a.shape)
return a.astype(np.float64)
class _Half:
"""Validated per-group integers and frequency rows of one half."""
def __init__(self, h, ng, name):
self.n = _exact_ints(h['n'], f'{name}.n')
self.s = _exact_ints(h['s'], f'{name}.s')
self.fg = np.asarray(h['fg']).astype(np.int64)
self.fd = np.asarray(h['fd']).astype(np.int64)
self.fc = _exact_ints(h['fc'], f'{name}.fc')
for arr, label in ((h['fg'], 'fg'), (h['fd'], 'fd')):
if np.asarray(arr).dtype.kind not in 'iu':
raise TypeError(f'{name}.{label} must be integer')
if len(self.n) != ng or len(self.s) != ng:
raise ValueError(f'{name}.n and {name}.s must have length n_groups={ng}')
if not (len(self.fg) == len(self.fd) == len(self.fc)):
raise ValueError(f'{name} frequency rows fg/fd/fc must have equal length')
if np.any(self.n < 0) or np.any(self.s < 0):
raise ValueError(f'{name}.n and {name}.s must be non-negative')
if len(self.fg):
if self.fg.min() < 0 or self.fg.max() >= ng or self.fd.min() < 1 or np.any(self.fc < 1):
raise ValueError(f'{name} rows need 0 <= group < n_groups, duration >= 1 and count >= 1')
ordered = (self.fg[1:] > self.fg[:-1]) | ((self.fg[1:] == self.fg[:-1]) & (self.fd[1:] > self.fd[:-1]))
if not np.all(ordered):
raise ValueError(f'{name} rows must be sorted by (group, duration) without duplicates')
ids = np.arange(ng)
self.start = np.searchsorted(self.fg, ids, side='left')
self.stop = np.searchsorted(self.fg, ids, side='right')
n_chk, s_chk = self._sums(ng)
if not (np.all(n_chk == self.n) and np.all(s_chk == self.s)):
raise ValueError(f'{name}: n/s are inconsistent with the frequency rows')
def _sums(self, ng):
fc = self.fc
small = fc.dtype != object and (not len(fc) or int(fc.max()) * int(self.fd.max()) * len(fc) < 2 ** 62)
if small:
n = np.zeros(ng, np.int64)
s = np.zeros(ng, np.int64)
np.add.at(n, self.fg, fc)
np.add.at(s, self.fg, fc * self.fd)
return n, s
n = [0] * ng
s = [0] * ng
for g, d, c in zip(self.fg.tolist(), self.fd.tolist(), [int(v) for v in fc.tolist()]):
n[g] += c
s[g] += c * d
return _obj(n), _obj(s)
def to_python_ints(self):
self.n, self.s, self.fc = _as_python_ints(self.n), _as_python_ints(self.s), _as_python_ints(self.fc)
return self
def rows(self, g):
return self.fd[self.start[g]:self.stop[g]], self.fc[self.start[g]:self.stop[g]]
[docs]
def joint_set(early, late):
"""Ascending ids of groups with at least one complete visit in both halves."""
nE = _exact_ints(early['n'], 'early.n')
nL = _exact_ints(late['n'], 'late.n')
return np.flatnonzero((nE >= 1) & (nL >= 1)).astype(np.int64)
# ------------------------------------------------------------------------------- production floats
def _sw_rows(values, masses, total, count):
"""Canonical S_w for each row of ``values`` (R x G), event masses over ``total`` and equal over ``count``.
Bitwise the PF ``_sw_pair`` recurrence: stable sort, exact-tie blocks, A recurrence and S sum
as sequential prefix sums along each row (adding +0.0 at non-block positions is exact).
"""
R, G = values.shape
if G <= 1:
return np.zeros(R), np.zeros(R)
order = np.argsort(values, axis=1, kind='stable')
sv = np.take_along_axis(values, order, axis=1)
sm = masses[order]
pos = np.arange(G)
start = np.ones((R, G), bool)
start[:, 1:] = sv[:, 1:] != sv[:, :-1]
end = np.ones((R, G), bool)
end[:, :-1] = start[:, 1:]
csum = np.cumsum(sm, axis=1)
before = csum - sm
next_end = np.minimum.accumulate(np.where(end, pos, G)[:, ::-1], axis=1)[:, ::-1]
block_mass = np.take_along_axis(csum, next_end, axis=1) - before
block_count = next_end - pos + 1
tv = float(total)
tc = float(count)
gap = np.zeros((R, G))
gap[:, 1:] = sv[:, 1:] - sv[:, :-1]
update = start.copy()
update[:, 0] = False
a_ev = np.cumsum(np.where(update, (_to_f64(before) / tv) * gap, 0.0), axis=1)
a_eq = np.cumsum(np.where(update, (pos / tc) * gap, 0.0), axis=1)
s_ev = np.cumsum(np.where(start, (_to_f64(block_mass) / tv) * a_ev, 0.0), axis=1)[:, -1]
s_eq = np.cumsum(np.where(start, (block_count / tc) * a_eq, 0.0), axis=1)[:, -1]
return 2.0 * s_ev, 2.0 * s_eq
def _canonical_sum(terms):
"""Sequential sum of the ascending-sorted terms (PF ``_canonical_sum``)."""
s = np.sort(np.asarray(terms, np.float64))
return float(np.cumsum(s)[-1]) if len(s) else 0.0
def _support_matrix(h, J):
"""Union support U (sorted durations) of the joint groups and cumulative counts C (|U| x G)."""
G = len(J)
idx = np.concatenate([np.arange(h.start[g], h.stop[g]) for g in J.tolist()])
ed = h.fd[idx]
ej = np.repeat(np.arange(G), (h.stop[J] - h.start[J]))
U = np.unique(ed)
cnt = np.zeros((len(U), G), dtype=h.fc.dtype)
cnt[np.searchsorted(U, ed), ej] = h.fc[idx]
return U, np.cumsum(cnt, axis=0)
def _bp_half(h, J, masses, total, chunk):
"""Integrate S_w(CDF) over the union support of one half (PF ``_bp_half``), both weightings."""
G = len(J)
U, C = _support_matrix(h, J)
support, maxd = len(U), int(U[-1])
if support < 2:
return 0.0, 0.0, support, maxd
dens = _to_f64(h.n[J])
s_ev = np.empty(support - 1)
s_eq = np.empty(support - 1)
step = max(1, int(chunk) // G)
for i0 in range(0, support - 1, step):
i1 = min(support - 1, i0 + step)
s_ev[i0:i1], s_eq[i0:i1] = _sw_rows(_to_f64(C[i0:i1]) / dens, masses, total, G)
gaps = (U[1:] - U[:-1]).astype(np.float64)
return float(np.cumsum(s_ev * gaps)[-1]), float(np.cumsum(s_eq * gaps)[-1]), support, maxd
def _w1_numerators(E, L, J):
"""Exact N_W1 = sum_i |C_E(u_i) n_L - C_L(u_i) n_E| (u_{i+1} - u_i) per joint group."""
out = []
for g in J.tolist():
dE, cE = E.rows(g)
dL, cL = L.rows(g)
u = np.union1d(dE, dL)
if len(u) < 2:
out.append(0)
continue
ce = np.zeros(len(u), dtype=cE.dtype)
cl = np.zeros(len(u), dtype=cL.dtype)
ce[np.searchsorted(u, dE)] = cE
cl[np.searchsorted(u, dL)] = cL
diff = np.abs(np.cumsum(ce) * L.n[g] - np.cumsum(cl) * E.n[g])
out.append(int(np.sum(diff[:-1] * np.diff(u))))
return _exact_ints(_obj(out), 'nw1')
def _same_distribution(h, g0, g):
"""Exact equality of two empirical distributions (same support, proportional counts)."""
a0, a1 = h.start[g0], h.stop[g0]
b0, b1 = h.start[g], h.stop[g]
if a1 - a0 != b1 - b0 or not np.array_equal(h.fd[a0:a1], h.fd[b0:b1]):
return False
return bool(np.all(h.fc[a0:a1] * h.n[g] == h.fc[b0:b1] * h.n[g0]))
[docs]
def status_code(G, b_zero, e_zero, b, e):
"""PF status: unavailable (empty joint set), indeterminate (floats contradict the exact identities
or are non-finite), finite (E > 0), undefined (B = E = 0) or divergent (B > 0, E = 0)."""
if G == 0:
return UNAVAILABLE
if not (np.isfinite(b) and np.isfinite(e)) or (b == 0.0) != bool(b_zero) or (e == 0.0) != bool(e_zero):
return INDETERMINATE
if not e_zero:
return FINITE
return UNDEFINED_00 if b_zero else DIVERGENT
def _combine(fo, io):
f = dict(zip(FLOAT_FIELDS, fo))
i = dict(zip(INT_FIELDS, (int(x) for x in io)))
out = {}
for cand, wt in COMBOS:
se = f[f'sep_{cand}_early_{wt}']
sl = f[f'sep_{cand}_late_{wt}']
b = 0.5 * (se + sl) if i['G_joint'] else np.nan
e = f[f'E_{cand}_{wt}']
b_zero = bool(i[f'{cand}_B_zero_early'] and i[f'{cand}_B_zero_late'])
e_zero = bool(i[f'{cand}_E_zero'])
code = status_code(i['G_joint'], b_zero, e_zero, b, e)
out[(cand, wt)] = {'B': b, 'E': e, 'ratio': b / e if code == FINITE else np.nan,
'separation_early': se, 'separation_late': sl,
'status_code': code, 'ratio_status': STATUS[code],
'B_exact_zero': b_zero if i['G_joint'] else None,
'E_exact_zero': e_zero if i['G_joint'] else None}
return out
def _check_guard_limit(guard_limit):
if isinstance(guard_limit, bool) or not isinstance(guard_limit, (int, np.integer)) or guard_limit > GUARD_LIMIT:
raise ValueError(f'guard_limit must be an integer <= GUARD_LIMIT = 2**62 (it may only lower the int64 '
f'bound), got {guard_limit!r}')
def _guard_ok(E, L, J, guard_limit):
"""Observed-bound guard computed in exact integers before any product."""
nE, nL = [int(v) for v in E.n[J].tolist()], [int(v) for v in L.n[J].tolist()]
sE, sL = [int(v) for v in E.s[J].tolist()], [int(v) for v in L.s[J].tolist()]
maxn = max(max(nE), max(nL))
maxs = max(max(sE), max(sL))
span = max(max(int(E.fd[E.stop[g] - 1]), int(L.fd[L.stop[g] - 1])) for g in J.tolist())
return not (maxn * maxs >= guard_limit or maxn * maxn * (span + 1) >= guard_limit
or maxn * maxn * 4 >= guard_limit)
[docs]
def separation(n_groups, early, late, *, guard_limit=GUARD_LIMIT, chunk_elements=1 << 20):
"""Production B/E for both candidates and weightings of one paired early/late dataset.
``early``/``late`` are mappings with ``n``, ``s`` (length ``n_groups``) and frequency rows
``fg, fd, fc`` sorted by (group, duration), e.g. from :func:`runs.half_stats`. Returns a dict
with ``path`` ('int64' or 'python_int'), ``joint_ids``, PF-ordered ``raw_float``/``raw_int``,
exact integer ``ntau``/``nw1``/``den`` and ``visit_mass`` (int64 arrays, or object arrays of
unclamped Python ints on the Python-integer path), per-group float ``dmu``/``w1``, the
``weights`` and per-group ``E_contributions`` that the canonical sums add, and per-combination
``combos[(candidate, weighting)]`` with B, E, ratio, half components and status.
``guard_limit`` may only lower :data:`GUARD_LIMIT` (``ValueError`` otherwise).
"""
_check_guard_limit(guard_limit)
ng = int(n_groups)
E = _Half(early, ng, 'early')
L = _Half(late, ng, 'late')
J = np.flatnonzero((E.n >= 1) & (L.n >= 1)).astype(np.int64)
G = len(J)
fo = np.full(12, np.nan)
io = np.zeros(13, np.int64)
io[0] = G
io[12] = 1
base = {'joint_ids': J, 'raw_float': fo, 'raw_int': io}
if G == 0:
z = np.zeros(0, np.int64)
zf = np.zeros(0)
base.update({'path': 'int64', 'visit_mass': z, 'visit_mass_total': 0, 'ntau': z, 'nw1': z, 'den': z,
'dmu': zf, 'w1': zf, 'mu_early': zf, 'mu_late': zf, 'weights': {'visit': zf, 'equal': zf},
'E_contributions': {c: zf for c in COMBOS}, 'combos': _combine(fo, io)})
return base
path = 'int64'
if not _guard_ok(E, L, J, guard_limit):
E.to_python_ints()
L.to_python_ints()
path = 'python_int'
nE, nL, sE, sL = E.n[J], L.n[J], E.s[J], L.s[J]
masses = nE + nL
total = int(np.sum(masses))
mu_e = _to_f64(sE) / _to_f64(nE)
mu_l = _to_f64(sL) / _to_f64(nL)
ntau = np.abs(sE * nL - sL * nE)
den = nE * nL
nw1 = _w1_numerators(E, L, J)
if path == 'python_int':
nw1 = _as_python_ints(nw1)
g0 = int(J[0])
io[1] = total if total <= _INT64_MAX else -1
io[2] = int(np.all(sE * nE[0] == sE[0] * nE))
io[3] = int(np.all(sL * nL[0] == sL[0] * nL))
io[4] = int(np.all(ntau == 0))
io[5] = int(all(_same_distribution(E, g0, g) for g in J.tolist()))
io[6] = int(all(_same_distribution(L, g0, g) for g in J.tolist()))
io[7] = int(np.all(nw1 == 0))
dmu = _to_f64(ntau) / _to_f64(den)
w1 = _to_f64(nw1) / _to_f64(den)
fo[0], fo[2] = (x[0] for x in _sw_rows(mu_e[None, :], masses, total, G))
fo[1], fo[3] = (x[0] for x in _sw_rows(mu_l[None, :], masses, total, G))
w_visit = _to_f64(masses) / float(total)
w_equal = 1.0 / float(G)
contrib = {('tau', 'visit'): w_visit * dmu, ('tau', 'equal'): w_equal * dmu,
('P', 'visit'): w_visit * w1, ('P', 'equal'): w_equal * w1}
fo[4] = _canonical_sum(contrib[('tau', 'visit')])
fo[5] = _canonical_sum(contrib[('tau', 'equal')])
fo[10] = _canonical_sum(contrib[('P', 'visit')])
fo[11] = _canonical_sum(contrib[('P', 'equal')])
fo[6], fo[8], io[8], io[10] = _bp_half(E, J, masses, total, chunk_elements)
fo[7], fo[9], io[9], io[11] = _bp_half(L, J, masses, total, chunk_elements)
base.update({'path': path, 'visit_mass': masses, 'visit_mass_total': total, 'ntau': ntau, 'nw1': nw1,
'den': den, 'dmu': dmu, 'w1': w1, 'mu_early': mu_e, 'mu_late': mu_l,
'weights': {'visit': w_visit, 'equal': np.full(G, w_equal)}, 'E_contributions': contrib,
'combos': _combine(fo, io)})
return base
# ------------------------------------------------------------------------------- Level 1 enclosures
def _dn(x):
y = math.nextafter(x, -math.inf)
return y if y > 0.0 else 0.0
def _up(x):
return math.nextafter(x, math.inf)
def _iv_int(x):
x = int(x)
f = float(x)
if -_TWO53 <= x <= _TWO53:
return f, f
return math.nextafter(f, -math.inf), math.nextafter(f, math.inf)
def _iv_div(num, den):
"""Enclosure of the exact non-negative rational num/den (den > 0)."""
if num == 0:
return 0.0, 0.0
a_lo, a_hi = _iv_int(num)
b_lo, b_hi = _iv_int(den)
return _dn(a_lo / b_hi), _up(a_hi / b_lo)
def _iv_add(alo, ahi, blo, bhi):
return _dn(alo + blo), _up(ahi + bhi)
def _iv_mul(alo, ahi, blo, bhi):
# operands enclose non-negative quantities with non-negative lower bounds
return _dn(alo * blo), _up(ahi * bhi)
def _dn_v(x):
y = np.nextafter(x, -np.inf)
return np.where(y > 0.0, y, 0.0)
def _up_v(x):
return np.nextafter(x, np.inf)
def _iv_int_v(x):
x = np.asarray(x)
f = _to_f64(x)
big = np.asarray((x > _TWO53) | (x < -_TWO53), bool)
return np.where(big, np.nextafter(f, -np.inf), f), np.where(big, np.nextafter(f, np.inf), f)
def _iv_div_v(num, den):
num = np.asarray(num)
a_lo, a_hi = _iv_int_v(num)
b_lo, b_hi = _iv_int_v(den)
zero = np.asarray(num == 0, bool)
return np.where(zero, 0.0, _dn_v(a_lo / b_hi)), np.where(zero, 0.0, _up_v(a_hi / b_lo))
def _iv_add_v(alo, ahi, blo, bhi):
return _dn_v(alo + blo), _up_v(ahi + bhi)
def _iv_mul_v(alo, ahi, blo, bhi):
return _dn_v(alo * blo), _up_v(ahi * bhi)
def _exact_order(num, den):
"""Per-row order of num[r, j]/den[j] that is exactly non-decreasing (ties in index order).
Fast path: PF's stable float-key sort; rows whose adjacent pairs violate the exact
cross-multiplied order are re-sorted with exact integer comparisons (PF aborted instead).
"""
R, G = num.shape
order = np.argsort(_to_f64(num) / _to_f64(den), axis=1, kind='stable')
sn = np.take_along_axis(num, order, axis=1)
sd = den[order]
bad = np.flatnonzero(np.any(np.asarray(sn[:, :-1] * sd[:, 1:] > sn[:, 1:] * sd[:, :-1], bool), axis=1))
for r in bad.tolist():
nr = [int(v) for v in num[r].tolist()]
dr = [int(v) for v in den.tolist()]
cmp = lambda x, y: (nr[x] * dr[y] > nr[y] * dr[x]) - (nr[x] * dr[y] < nr[y] * dr[x])
order[r] = sorted(range(G), key=cmp_to_key(cmp))
return order, len(bad)
def _sw_iv_rows(num, den, masses, total, count):
"""Level-1 interval S_w for each row of exact rationals num[r, j]/den[j] (port of PF ``_sw_iv``)."""
R, G = num.shape
z = np.zeros(R)
if G <= 1:
return z, z.copy(), z.copy(), z.copy(), 0
order, reordered = _exact_order(num, den)
sn = np.take_along_axis(num, order, axis=1)
sd = den[order]
sm = masses[order]
pos = np.arange(G)
start = np.ones((R, G), bool)
start[:, 1:] = np.asarray(sn[:, 1:] * sd[:, :-1] != sn[:, :-1] * sd[:, 1:], bool)
end = np.ones((R, G), bool)
end[:, :-1] = start[:, 1:]
csum = np.cumsum(sm, axis=1)
before = csum - sm
next_end = np.minimum.accumulate(np.where(end, pos, G)[:, ::-1], axis=1)[:, ::-1]
block_mass = np.take_along_axis(csum, next_end, axis=1) - before
block_count = next_end - pos + 1
a_ev = [np.zeros(R), np.zeros(R)]
a_eq = [np.zeros(R), np.zeros(R)]
s_ev = [np.zeros(R), np.zeros(R)]
s_eq = [np.zeros(R), np.zeros(R)]
prev = np.full(R, -1)
for k in range(G):
rows = np.flatnonzero(start[:, k])
if not len(rows):
continue
if k > 0:
p = prev[rows]
glo, ghi = _iv_div_v(sn[rows, k] * sd[rows, p] - sn[rows, p] * sd[rows, k], sd[rows, k] * sd[rows, p])
for acc, wnum, wden in ((a_ev, before[rows, k], total), (a_eq, np.full(len(rows), k), count)):
wlo, whi = _iv_div_v(wnum, wden)
plo, phi = _iv_mul_v(wlo, whi, glo, ghi)
acc[0][rows], acc[1][rows] = _iv_add_v(acc[0][rows], acc[1][rows], plo, phi)
for acc, a, wnum, wden in ((s_ev, a_ev, block_mass[rows, k], total), (s_eq, a_eq, block_count[rows, k], count)):
wlo, whi = _iv_div_v(wnum, wden)
plo, phi = _iv_mul_v(wlo, whi, a[0][rows], a[1][rows])
acc[0][rows], acc[1][rows] = _iv_add_v(acc[0][rows], acc[1][rows], plo, phi)
prev[rows] = k
return 2.0 * s_ev[0], 2.0 * s_ev[1], 2.0 * s_eq[0], 2.0 * s_eq[1], reordered
[docs]
def level1(n_groups, early, late, result, *, chunk_elements=1 << 20):
"""Level-1 enclosures of the exact B and E of every combination (port of PF ``exact_verify.level1``).
Uses the production ``result`` for the joint set, the integer path and the exact numerators.
Inside the int64 guard all integer work is int64; on ``result['path'] == 'python_int'`` it
uses Python integers with the exact, unclamped ``ntau/nw1/den``. Returns ``None`` for an
empty joint set, else ``{'path', 'enclosures': {combo: (B_lo, B_hi, E_lo, E_hi)}, 'halves',
'reordered_rows'}`` where ``reordered_rows`` counts float-sort orders repaired exactly.
"""
J = np.asarray(result['joint_ids'], np.int64)
G = len(J)
if G == 0:
return None
ng = int(n_groups)
E = _Half(early, ng, 'early')
L = _Half(late, ng, 'late')
python_int = result['path'] == 'python_int'
if python_int:
E.to_python_ints()
L.to_python_ints()
masses = E.n[J] + L.n[J]
total = int(np.sum(masses))
halves = {}
reordered = 0
for tag, h in (('early', E), ('late', L)):
lo_ev, hi_ev, lo_eq, hi_eq, r = _sw_iv_rows(h.s[J][None, :], h.n[J], masses, total, G)
reordered += r
halves[('tau', 'visit', tag)] = (float(lo_ev[0]), float(hi_ev[0]))
halves[('tau', 'equal', tag)] = (float(lo_eq[0]), float(hi_eq[0]))
for tag, h in (('early', E), ('late', L)):
U, C = _support_matrix(h, J)
acc = {'visit': (0.0, 0.0), 'equal': (0.0, 0.0)}
step = max(1, int(chunk_elements) // G)
for i0 in range(0, len(U) - 1, step):
i1 = min(len(U) - 1, i0 + step)
lo_ev, hi_ev, lo_eq, hi_eq, r = _sw_iv_rows(C[i0:i1], h.n[J], masses, total, G)
reordered += r
for i in range(i1 - i0):
glo, ghi = _iv_int(int(U[i0 + i + 1]) - int(U[i0 + i]))
for wt, slo, shi in (('visit', lo_ev[i], hi_ev[i]), ('equal', lo_eq[i], hi_eq[i])):
plo, phi = _iv_mul(float(slo), float(shi), glo, ghi)
acc[wt] = _iv_add(acc[wt][0], acc[wt][1], plo, phi)
halves[('P', 'visit', tag)] = acc['visit']
halves[('P', 'equal', tag)] = acc['equal']
ntau, nw1, den = result['ntau'], result['nw1'], result['den']
tv = te = pv = pe = (0.0, 0.0)
q = _iv_div(1, G)
for j in range(G):
w = _iv_div(int(masses[j]), total)
d = _iv_div(int(ntau[j]), int(den[j]))
tv = _iv_add(*tv, *_iv_mul(*w, *d))
te = _iv_add(*te, *_iv_mul(*q, *d))
d = _iv_div(int(nw1[j]), int(den[j]))
pv = _iv_add(*pv, *_iv_mul(*w, *d))
pe = _iv_add(*pe, *_iv_mul(*q, *d))
denominators = {('tau', 'visit'): tv, ('tau', 'equal'): te, ('P', 'visit'): pv, ('P', 'equal'): pe}
enclosures = {}
for cand, wt in COMBOS:
lo, hi = _iv_add(*halves[(cand, wt, 'early')], *halves[(cand, wt, 'late')])
enclosures[(cand, wt)] = (0.5 * lo, 0.5 * hi) + denominators[(cand, wt)]
return {'path': result['path'], 'enclosures': enclosures, 'halves': halves, 'reordered_rows': int(reordered)}
# ------------------------------------------------------------------------------- Level 2b exact
[docs]
def exact_components(n_groups, early, late, *, max_seconds=None, max_work=None):
"""Exact ``Fraction`` B and E of every combination (port of PF ``exact_components('fraction')``).
``max_work`` bounds the predicted number of exact block updates (joint groups x support
points); ``max_seconds`` bounds the wall time of the whole evaluation, checked throughout and
once more before returning. Exceeding either raises :class:`BudgetExhausted`.
Returns ``None`` for an empty joint set, else ``{combo: (B, E)}``. With neither bound the
evaluation is unbounded: a low-level compatibility default. The public assessment route
(:mod:`.assessment`) only calls it under an explicit finite budget.
"""
deadline = None if max_seconds is None else time.monotonic() + float(max_seconds)
def tick():
if deadline is not None and time.monotonic() > deadline:
raise BudgetExhausted(f'exact evaluation exceeded max_seconds={max_seconds}')
ng = int(n_groups)
E = _Half(early, ng, 'early').to_python_ints()
L = _Half(late, ng, 'late').to_python_ints()
nE = [int(x) for x in E.n.tolist()]
sE = [int(x) for x in E.s.tolist()]
nL = [int(x) for x in L.n.tolist()]
sL = [int(x) for x in L.s.tolist()]
J = [g for g in range(ng) if nE[g] >= 1 and nL[g] >= 1]
if not J:
return None
G = len(J)
rows = {}
for tag, h in (('E', E), ('L', L)):
rows[tag] = {g: list(zip(h.rows(g)[0].tolist(), [int(c) for c in h.rows(g)[1].tolist()])) for g in J}
if max_work is not None:
support = sum(len({d for g in J for d, _ in rows[t][g]}) for t in ('E', 'L'))
if G * (support + 2) > max_work:
raise BudgetExhausted(f'exact work {G * (support + 2)} > max_work {max_work}')
tick()
masses = [nE[g] + nL[g] for g in J]
T = sum(masses)
def sw(pairs):
order = sorted(range(G), key=lambda j: Fraction(pairs[j][0], pairs[j][1]))
a_ev = a_eq = s_ev = s_eq = Fraction(0)
cum_ev = cum_eq = 0
prev = None
j = 0
while j < G:
tick()
a0 = order[j]
m_ev = m_eq = 0
k = j
while k < G and pairs[order[k]][0] * pairs[a0][1] == pairs[a0][0] * pairs[order[k]][1]:
m_ev += masses[order[k]]
m_eq += 1
k += 1
if prev is not None:
gap = Fraction(pairs[a0][0] * pairs[prev][1] - pairs[prev][0] * pairs[a0][1], pairs[a0][1] * pairs[prev][1])
a_ev += Fraction(cum_ev, T) * gap
a_eq += Fraction(cum_eq, G) * gap
s_ev += Fraction(m_ev, T) * a_ev
s_eq += Fraction(m_eq, G) * a_eq
cum_ev += m_ev
cum_eq += m_eq
prev = a0
j = k
return 2 * s_ev, 2 * s_eq
t_e = sw([(sE[g], nE[g]) for g in J])
t_l = sw([(sL[g], nL[g]) for g in J])
sep = {}
for tag, n in (('E', nE), ('L', nL)):
events = sorted((d, j, c) for j, g in enumerate(J) for d, c in rows[tag][g])
cum = [0] * G
acc_ev = acc_eq = Fraction(0)
i = 0
while i < len(events):
tick()
d = events[i][0]
k = i
while k < len(events) and events[k][0] == d:
cum[events[k][1]] += events[k][2]
k += 1
if k < len(events):
gap = events[k][0] - d
s_ev, s_eq = sw([(cum[j], n[J[j]]) for j in range(G)])
acc_ev += s_ev * gap
acc_eq += s_eq * gap
i = k
sep[tag] = (acc_ev, acc_eq)
et_ev = et_eq = ep_ev = ep_eq = Fraction(0)
for j, g in enumerate(J):
tick()
dnum = abs(sE[g] * nL[g] - sL[g] * nE[g])
den = nE[g] * nL[g]
ce, cl = dict(rows['E'][g]), dict(rows['L'][g])
support = sorted(set(ce) | set(cl))
cum_e = cum_l = w1 = 0
for i, u in enumerate(support):
cum_e += ce.get(u, 0)
cum_l += cl.get(u, 0)
if i + 1 < len(support):
w1 += abs(cum_e * nL[g] - cum_l * nE[g]) * (support[i + 1] - u)
et_ev += Fraction(masses[j], T) * Fraction(dnum, den)
et_eq += Fraction(1, G) * Fraction(dnum, den)
ep_ev += Fraction(masses[j], T) * Fraction(w1, den)
ep_eq += Fraction(1, G) * Fraction(w1, den)
half = Fraction(1, 2)
tick()
return {('tau', 'visit'): ((t_e[0] + t_l[0]) * half, et_ev), ('tau', 'equal'): ((t_e[1] + t_l[1]) * half, et_eq),
('P', 'visit'): ((sep['E'][0] + sep['L'][0]) * half, ep_ev),
('P', 'equal'): ((sep['E'][1] + sep['L'][1]) * half, ep_eq)}
[docs]
def reliability(result, enclosures=None, exact=None):
"""ARITHMETIC reliability fields per combination (rounding error of the exact value only).
``enclosures``: :func:`level1` output (or ``None``); ``exact``: :func:`exact_components`
output (or ``None``). Relative errors are NaN with reason ``'exact_zero'`` when the exact
value is zero (absolute error and enclosure widths remain), or ``'exact_not_computed'``.
"""
out = {}
for combo in COMBOS:
c = result['combos'][combo]
rec = {'reliability_kind': RELIABILITY_KIND}
box = enclosures['enclosures'][combo] if enclosures else (np.nan,) * 4
rec.update({'B_enclosure_lo': box[0], 'B_enclosure_hi': box[1], 'E_enclosure_lo': box[2],
'E_enclosure_hi': box[3], 'B_enclosure_width': box[1] - box[0], 'E_enclosure_width': box[3] - box[2],
# False is not an error: the canonical production float may round just outside the
# enclosure of the exact value (benign cancellation); fronts use the enclosure.
'production_inside_enclosure': bool(enclosures is not None and box[0] <= c['B'] <= box[1]
and box[2] <= c['E'] <= box[3])})
for q, k in (('B', 0), ('E', 1)):
if exact is None or not np.isfinite(c[q]):
rec.update({f'{q}_abs_error_vs_exact': np.nan, f'{q}_rel_error_vs_exact': np.nan,
f'{q}_rel_error_reason': 'exact_not_computed'})
continue
x = exact[combo][k]
err = abs(Fraction(float(c[q])) - x)
rec[f'{q}_abs_error_vs_exact'] = float(err)
rec[f'{q}_rel_error_vs_exact'] = float(err / x) if x != 0 else np.nan
rec[f'{q}_rel_error_reason'] = '' if x != 0 else 'exact_zero'
out[combo] = rec
return out
COVERAGE = ('occupied_groups_early', 'occupied_groups_late', 'complete_groups_early', 'complete_groups_late',
'joint_groups', 'complete_visits_early', 'complete_visits_late', 'joint_complete_visits_early',
'joint_complete_visits_late', 'complete_occupancy_early', 'complete_occupancy_late',
'joint_complete_occupancy_early', 'joint_complete_occupancy_late', 'accepted_occupancy_early',
'accepted_occupancy_late', 'joint_accepted_occupancy_early', 'joint_accepted_occupancy_late',
'complete_entities_early', 'complete_entities_late', 'joint_complete_entities_early',
'joint_complete_entities_late', 'early_only_groups', 'late_only_groups', 'early_only_complete_visits',
'late_only_complete_visits', 'early_only_complete_occupancy', 'late_only_complete_occupancy',
'early_only_accepted_occupancy', 'late_only_accepted_occupancy', 'occupied_zero_complete_groups_early',
'occupied_zero_complete_groups_late', 'one_frame_visits_early', 'one_frame_visits_late',
'joint_one_frame_visits_early', 'joint_one_frame_visits_late', 'left_only_early', 'left_only_late',
'right_only_early', 'right_only_late', 'both_censored_early', 'both_censored_late',
'rejected_occupancy_early', 'rejected_occupancy_late', 'rejected_runs_early', 'rejected_runs_late',
'frames_changed_early', 'frames_changed_late', 'accepted_to_accepted_changes_early',
'accepted_to_accepted_changes_late', 'rejected_to_accepted_changes_early',
'rejected_to_accepted_changes_late', 'overlapping_candidate_pairs_early', 'overlapping_candidate_pairs_late')
[docs]
def coverage(early, late, joint_ids):
"""Exact integer coverage/censoring/bridge accounting of one paired dataset (PF ``coverage``).
Bridge-change fields are ``None`` when a half's counts are unknown (see :func:`runs.half_stats`,
whose ``bridge_changes_reason`` says why); they are never reported as 0.
"""
J = np.zeros(len(early['n']), bool)
J[np.asarray(joint_ids, np.int64)] = True
out = {}
for tag, h, other in (('early', early, late), ('late', late, early)):
occ_g = h['occ'] > 0
comp_g = h['n'] > 0
only = comp_g & (other['n'] == 0)
out[f'occupied_groups_{tag}'] = int(occ_g.sum())
out[f'complete_groups_{tag}'] = int(comp_g.sum())
out[f'complete_visits_{tag}'] = int(h['n'].sum())
out[f'joint_complete_visits_{tag}'] = int(h['n'][J].sum())
out[f'complete_occupancy_{tag}'] = int(h['s'].sum())
out[f'joint_complete_occupancy_{tag}'] = int(h['s'][J].sum())
out[f'accepted_occupancy_{tag}'] = int(h['occ'].sum())
out[f'joint_accepted_occupancy_{tag}'] = int(h['occ'][J].sum())
out[f'complete_entities_{tag}'] = int(h['ogrp'].any(axis=1).sum())
out[f'joint_complete_entities_{tag}'] = int(h['ogrp'][:, J].any(axis=1).sum())
out[f'{tag}_only_groups'] = int(only.sum())
out[f'{tag}_only_complete_visits'] = int(h['n'][only].sum())
out[f'{tag}_only_complete_occupancy'] = int(h['s'][only].sum())
out[f'{tag}_only_accepted_occupancy'] = int(h['occ'][only].sum())
out[f'occupied_zero_complete_groups_{tag}'] = int((occ_g & ~comp_g).sum())
out[f'one_frame_visits_{tag}'] = int(h['one'].sum())
out[f'joint_one_frame_visits_{tag}'] = int(h['one'][J].sum())
out[f'left_only_{tag}'] = int(h['left_only'].sum())
out[f'right_only_{tag}'] = int(h['right_only'].sum())
out[f'both_censored_{tag}'] = int(h['both'].sum())
out[f'rejected_occupancy_{tag}'] = h['rejected_occupancy_frames']
out[f'rejected_runs_{tag}'] = h['rejected_runs']
for key in BRIDGE_FIELDS:
out[f'{key}_{tag}'] = h[key]
out['joint_groups'] = int(J.sum())
return {key: out[key] for key in COVERAGE}
def _same_ids(a, b):
a, b = np.asarray(a), np.asarray(b)
return a.shape == b.shape and all(x == y for x, y in zip(a.tolist(), b.tolist()))
[docs]
def check_halves(early_runs, late_runs):
"""Refuse (``ValueError``) an early/late pair that is not the two halves of the same record(s).
Both tables must share preprocessing, and source by source (in order) record identity, the
declared source population (both undeclared, or equal), entities, ``frame_interval`` and
``time_unit``; the early window must end at or before the
late start (not swapped/overlapping), and the two windows must be ``(0, s)`` and
``(s, n_frames)`` of that record, i.e. :meth:`record.LabelRecord.half_windows`. Duration
distributions from different sampling or different records are never paired.
"""
if early_runs.preprocessing != late_runs.preprocessing:
raise ValueError(f'early/late preprocessing differ: {early_runs.preprocessing!r} vs {late_runs.preprocessing!r}')
es, ls = early_runs.sources, late_runs.sources
if len(es) != len(ls) or not _same_ids(early_runs.entity_ids, late_runs.entity_ids):
raise ValueError('early/late record identity differs: different sources or entity identities')
for a, b in zip(es, ls):
if (a.record_id, a.entity_offset, a.n_entities) != (b.record_id, b.entity_offset, b.n_entities):
raise ValueError(f'early/late record identity differs: {a.record_id!r} vs {b.record_id!r}')
if getattr(a, 'source', None) != getattr(b, 'source', None):
raise ValueError(f'{a.record_id!r}: early/late declared source populations differ ({a.source!r} vs '
f'{b.source!r}); the two halves must come from one declared source')
if a.frame_interval != b.frame_interval:
raise ValueError(f'{a.record_id!r}: early/late frame_interval differ ({a.frame_interval!r} vs '
f'{b.frame_interval!r}); never pool durations from incompatible sampling')
if a.time_unit != b.time_unit:
raise ValueError(f'{a.record_id!r}: early/late time_unit differ ({a.time_unit!r} vs {b.time_unit!r})')
(e0, e1), (l0, l1) = ((int(x) for x in a.window), (int(x) for x in b.window))
if e1 > l0:
raise ValueError(f'{a.record_id!r}: early window {(e0, e1)} does not end at or before the late window '
f'{(l0, l1)} (swapped or overlapping)')
n = a.record_frames
if n is None or n != b.record_frames:
raise ValueError(f'{a.record_id!r}: record length unknown or different ({n!r} vs {b.record_frames!r}); '
'cannot verify that the windows are halves of one record')
if (e0, e1, l1) != (0, l0, n):
raise ValueError(f'{a.record_id!r}: windows {(e0, e1)} and {(l0, l1)} are not the two halves (0, s), '
f'(s, {n}) of the record window (0, {n})')
if a.record_fingerprint is None or b.record_fingerprint is None:
raise ValueError(f'{a.record_id!r}: record content identity unknown; cannot verify that both halves '
'come from one record')
if a.record_fingerprint != b.record_fingerprint:
raise ValueError(f'early/late record identity differs: {a.record_id!r} names two records with '
'different content (labels, identities, segments or sampling)')
[docs]
def assess(early_runs, late_runs, n_groups, *, guard_limit=GUARD_LIMIT, with_level1=True):
"""Convenience: half statistics, production B/E, Level-1 enclosures, coverage and reliability.
The pair is validated first by :func:`check_halves`, and ``guard_limit`` may only lower
:data:`GUARD_LIMIT`; both refuse with ``ValueError``.
"""
check_halves(early_runs, late_runs)
_check_guard_limit(guard_limit)
early = half_stats(early_runs, n_groups)
late = half_stats(late_runs, n_groups)
result = separation(n_groups, early, late, guard_limit=guard_limit)
l1 = level1(n_groups, early, late, result) if with_level1 and len(result['joint_ids']) else None
return {'early': early, 'late': late, 'result': result, 'level1': l1,
'coverage': coverage(early, late, result['joint_ids']), 'reliability': reliability(result, l1)}