residual.py
Cited on: PRINCIPAL.
To run it: save residual.py, with module_law.py, origin.py, principal.py beside it; needs numpy, scipy.
"""residual.py -- the real fields' remainder after the term at s = 0,
read one class character at a time.
QUESTION. origin.py reads, over a population of quadratic fields with
h+ >= 2, the level of a class C counted in prime powers, PI_C/(PI/h+),
against E_C, the explicit formula's term at s = 0. What is left on a
cell C is, per field, an identity of orthogonality:
level_C - 1 - E_C = sum_{chi != 1} conj chi(C) r_chi / PI(x, 1),
r_chi = PI(x, chi) - [-m (loglog x + gamma) + log L*(0, chi)],
the nontrivial zeros' sum of chi up to O(1/(x log x)). Over real fields
it is small at h+ = 2 and large at h+ >= 4. Is that growth the count of
nontrivial characters, each kind of character carrying a fixed mean
rho = r_chi/PI(x, 1), or does a character's own mean move with h+ or
with |D|?
KINDS. Real even (chi(k0) = +1, or no k0), real odd (chi(k0) = -1),
imaginary; each split genus (real-valued) or complex, the split printed
and not fitted.
HAND, from origin.py's F2 at 250. h+ = 2 with k0: one odd character,
rho_o = T = -K0, about +0.008; h+ = 3 (no k0, two complex even
characters): rho_e about -0.012 +- 0.007; h+ = 4 elementary with k0,
T = rho_e + 2 rho_o, K0 = rho_e - 2 rho_o: rho_e about +0.009, rho_o
about +0.019. So rho_o at h+ >= 4, at larger |D|, is about twice
h+ = 2's: a fixed mean pooled over |D| should fail, and |D| is where it
moves.
PREDICTIONS, frozen before the engine.
R1 printed: rho by kind x (genus, complex) x |D| band x h band, at 250
and 10^4; and the fields by h band x |D| band.
R2 printed: the field residual on T regressed on (h+ - 1)/PI(x, 1) at
every cut, real and imaginary.
R3 predicted: rho_o rises with log|D| at 250, slope > 0 by 3 se.
KILL FIXED. "A fixed mean per character of each kind carries the field
residual's dependence on h+": the kind means fitted over every real
field, a field's prediction the sum over its characters of conj
chi(C) times its kind's mean; dies if at 250 the T or K0 residual of
any h band (2, 3, 4, >= 5, and h >= 4 elementary 2 or not) departs
from its prediction by more than 3 of its own se. KILL FIXED|D|: the
same with the kind means fitted within each |D| band.
C1w every genus character of every field to |D| = C1_D, D = d1 d2,
signed: L*(0) real and positive, as the closed form to 1e-8.
C5 complex characters, the known answer: Lambda(2) from the kernels
W(y; 2) = int_1^oo phi(yt) (t + t^-2) dt, phi = 4 K0(2t) (even),
4t K0(2t) (odd), e^-t (imaginary), against L(2, chi)'s Euler product
over p < 10^4, every nontrivial character to |D| = C1_D, relative
gap below 2.5/(P ln P). Its break: the Euler product with
chi(P bar) = chi(P) fails at the complex characters (> 1e-3) and
passes at the genus ones.
C6 the identity: the characters' sum reproduces origin.py's cell
residual per field to 1e-10.
C7 THE KILL'S CONTROL: each rho replaced by its kind's mean plus noise
of its kind's spread, KILL FIXED misses; by the kind's mean times
h+/4, it fires.
THE KERNELS AT s = 2. With the sum Lambda~(s) = sum_n b_n int_0^oo
phi(nt/A) t^s dt/t = A^s G(s) L(s, chi), G = Gamma(s/2)^2,
Gamma((s+1)/2)^2, Gamma(s), and Theta(1/t) = t Theta(t),
Lambda~(s) = sum_n b_n W(n/A; s), W(y; s) = int_1^oo phi(yt)(t^(s-1) +
t^-s) dt, which at s = 1 is origin.py's W. At s = 2, z = 2y:
W_e = 4 K1(z)/z + 4z int_z^oo K0(v) v^-2 dv,
W_o = 2 I2(z)/z^2 + 2z int_z^oo K0(v)/v dv,
I2 = int_z^oo v^2 K0 = z^2 K1(z) + z K0(z) + int_z^oo K0,
W_c = (1 + y) e^-y / y^2 + e^-y - y E1(y),
and L(2) = sum/A^2 (even, imaginary), sum/(A^2 pi/4) (odd: Gamma(3/2)^2
= pi/4).
CHANGED AFTER THE FIRST FULL RUN, labelled so. C5's break failed its
frozen bar (complex least gap 6.53e-05 against > 1e-3, a hand guess at
the break's size), yet sits above C5's own bar at every complex
character, so C5 catches the break wherever it can act: the break's bar
is now C5's. Its genus half cannot fail, the break being the identity
on real characters. The pooled kind means, computed, became a print.
ADDED AT THE AUDIT, a check (C8): L*(0, chi) real and positive at every
nontrivial character of both populations, complex ones included (C1w
reads genus characters only); the trivial row is a positive theta sum
by construction, not zeta_K's L*(0), and is left out.
ADDED AFTER THE FIRST READING, frozen before its code (section C, real
fields, cut 250): the step in rho_o from h+ = 2 to h+ >= 4.
S0 printed: the step's targets, rho_o over the odd characters at h = 2
and at h >= 4; each kind's sd per h band.
C9u every complex character's conjugate is in its field, every genus
character stands alone: a conjugate pair is one unit, rho_chi-bar =
conj rho_chi. C9a the seeded kill reproduces KILL FIXED on the data.
C9 THE STEP'S CONTROL, 200 seeds, one draw per unit: plant FIXED (each
kind its pooled mean plus noise of its sd) and STEP (odd at the
step's targets, even at its pooled mean). STEP's z median above 1.5
under STEP, |median| below 0.5 under FIXED, FIXED's STEP sd within
25% of STEP's se, KILL FIXED firing below 0.10 under FIXED. Predicted:
KILL FIXED fires at 0.10 to 0.45 of the STEP seeds.
Z2 printed: the sum of z^2 over the T and K0 cells of h = 2, 3, 4, >= 5
(seven), and over the four T cells; p against the FIXED seeds, power
at their 95th.
KILL STEP. STEP = mean rho_o over units at h >= 4 minus at h = 2,
against the FIXED seeds' sd: above 3 STANDS, below 1 KILLED, between
UNDECIDED; if UNDECIDED, the real population widened to |D| <= 48000
(--dreal 48000), the same kill.
Not in the freeze, written with the code: C9u, C9a, and C9's clause
that KILL FIXED fires below 0.10 under FIXED.
CHANGED AFTER THE WIDENED RUN, labelled so. The STEP plant took the
population's own targets, so at 48000 it planted the widened step and
C9 failed (median 0.89); it now plants NAMED, the targets at 16000, at
either population. The prediction 0.10 to 0.45 is the 16000 one. The
line "the named step: observed ... sd from it" was added with it. NAMED
is read from the 16000 fields, which the 48000 population contains, so
its distance at 48000 is in sample, a bias toward the step.
ADDED AT THE AUDIT: R3u (R3 over units, a conjugate pair once, R3 having
counted each pair twice); each |D| band's STEP printed as its distance
from the named step; a band with fewer than two odd characters a side
skipped.
THE PARITY IDENTITY, by hand, not yet read. Over a real field with k0,
the even characters sum chi(C) to (h+/2)[C in {1, k0}], the odd ones to
(h+/2)([C = 1] - [C = k0]), so with PI_C the count in C:
odd mean of PI(x, chi) = PI_1 - PI_k0,
even nontrivial mean of PI(x, chi) = ((h+/2)(PI_1 + PI_k0) - PI)/(h+/2 - 1).
The even line holds at h+ >= 4, there being no even nontrivial
character at h+ = 2. The odd mean of rho is then (PI_1 - PI_k0 - the
odd mean of log L*(0, chi))/PI(x, 1): a fixed rho_o asks that
difference of absolute counts, the small primes the trivial class holds
over k0, not to scale with h+. And by the functional equation
log L*(0, chi) is (1/2) log D + log L(1, chi) plus a constant set by the
parity, so every rho carries -(1/2) log D/PI(x, 1) against whatever the
counts do with |D|.
FINDINGS, the real fields at the default |D| <= 48000 (12599 fields,
63511 nontrivial characters, 47189 units), the imaginary at |D| <= 4000
(1208 fields, 22385 characters).
F0 CONTROLS. C1w: 1177 genus characters to |D| <= 1500, L*(0) real
and positive at every one, the closed form to 9.18e-14. C5: L(2,
chi) from the s = 2 kernels against the Euler product, genus 1177
worst 4.43e-07, complex 4632 worst 4.13e-07, bar 2.72e-05; the
break moves every complex character by 6.53e-05 or more. C6 to
6.33e-15 (724 real fields), 3.50e-14 (902 imaginary). C7: the
fixed plant misses (worst +2.0 se), the h+/4 plant fires (-28.6
se). C8: L*(0) real and positive at all 85896 nontrivial
characters.
F1 KILL FIXED MISSED, real: pooled rho_e +0.0133, rho_o +0.0091 at
250 put every h band within 3 se, worst -2.8 se (h >= 4 elementary
2, K0 -0.0118 against -0.0025; its T +2.6 se); h = 2 T +0.0077
against +0.0101, h >= 5 T +0.1247 against +0.1245. KILL FIXED|D|
MISSED, worst +2.8 se (h >= 4 elementary 2, T). A fixed mean per
parity is not killed at its bar, and carries most of the
residual's growth with h+; the kill's control never planted the
step the hand reading named (rho_o from h = 2 to h >= 4, flat past
4): F6 prices it.
F2 R3 MISSED, the other way: rho_o on log|D| at 250, slope -0.0035
+- 0.0007 (-5.1 se, 32585 odd characters taken as independent);
R3u, over units, -0.0030 +- 0.0008 (-3.9 se, 25484). The kind
means by |D| band at 250, odd +0.0151, +0.0162, +0.0083, even
-0.0084, -0.0016, +0.0153: |D| moves both, the odd falling by
its slope, the even rising band on band, and pooled the even
mean sits above the odd. At 10^4 every
real R1 cell lies within 0.005 of 0.
F3 R2, the mean of r_chi in counts: real slope +0.63 +- 0.02 at 250,
+1.14, +1.30, +1.62, +1.67 at 400 to 2500, +1.86 +- 0.09 at 10^4,
intercepts +0.0043 +- 0.0026, +0.0099 +- 0.0022, +0.0082 +- 0.0018
at 250, 400, 630, within 0.003 of 0 from 1000 on; imaginary +0.88
+- 0.04 at 250 with intercept -0.172 +- 0.013, slope -0.16 +- 0.19
and intercept +0.001 at 10^4.
F4 printed: the imaginary arm of KILL FIXED fires at h = 4 (T -0.0669
+- 0.0237 against +0.0144, -3.4 se, 54 fields) and at h >= 4
elementary 2 (-3.3 se, 37 fields), the other bands within 0.5 se.
F5 the confound, counted: 3793 of the 4271 real fields with h+ >= 5,
and 3338 of the 4145 with h+ = 2, have |D| > 8000.
F6 THE STEP: targets +0.0075 (3219 odd characters at h = 2), +0.0093
(29366 at h >= 4). KILL STEP UNDECIDED at its frozen bar: +0.0023
+- 0.0020 over units (h >= 4 +0.0098, 22265 units; h = 2 +0.0075,
3219), +1.15 against the FIXED seeds' sd 0.0020, and 2.54 sd below
the named step +0.0073, which fires KILL FIXED at 0.735 of its
seeds and passes Z2's 95th at 0.915 (the print's "MISSED (0.735)"
is the prediction frozen at 16000). Z2 reads 7.12, p 0.215; Z2_T
5.59, p 0.135. NAMED is read from the |D| <= 16000 fields, which
this population contains: in sample, a bias toward the step. By
|D| band the step reads +3.2, -0.3, +0.8 se, from the named step
+1.8, -1.5, -2.5 se: past 8000, +0.0017 +- 0.0023 (2584 odd units
at h = 2), the band sits furthest below it, where at 16000 the
band 8000 to 16000 read +0.0070 +- 0.0048 (537), on it; the
shortfall is the fields past 16000's.
TIERS. F1 to F6 are observations over these populations.
RUN RECORD. A rehearsal at --d 600: 6/7 (C7's h+/4 plant at 134 fields
lacks power, -1.9 se). The first full run 36.7 s, peak commit 131.9 MB,
6/7 (C5's break as frozen); after the change 7/7, 36.7 s, 132.3 MB,
under a memory guard (512 MB). At the audit, with C8: 8/8, 36.9 s,
132.0 MB; C8 then narrowed to the nontrivial characters: 8/8, 37.1 s,
132.2 MB. With section C: a rehearsal at --d 600 9/11 (its cells too
small for C7 and C9); 11/11, 38.5 s, 135.2 MB; --dreal 48000 10/11, 114.2
s, 189.0 MB (C9 as frozen). After the change: 11/11, 38.2 s, 134.6 MB;
--dreal 48000 11/11, 113.4 s, 189.0 MB. At the audit (R3u, the band
distances): 11/11, 38.3 s, 135.4 MB; --dreal 48000 11/11, 114.1 s,
190.5 MB. The real default then moved from 16000 to 48000: 11/11, 112.5
s, 190.4 MB.
The reading at |D| <= 16000 (--dreal 16000; 4107 real fields, 16640
characters): KILL FIXED missed, pooled rho_e +0.0025, rho_o +0.0146,
worst +2.3 se (h >= 4 elementary 2, T); R3 +0.0009 +- 0.0012 (8751
characters); R2 real +0.42 +- 0.04 at 250, +1.43 +- 0.18 at 10^4. The
step's targets +0.0082 (1172), +0.0155 (7579), named from here; KILL
FIXED fired at 0.175 of the STEP seeds (predicted 0.10 to 0.45: HELD),
0.005 of the FIXED; KILL STEP UNDECIDED, +0.0073 +- 0.0032, +2.38
against the FIXED sd 0.0031; Z2 15.21, p 0.015 (power 0.53); by |D|
band +3.2, -0.3, +1.5 se; R3u +0.9 se.
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import sys
import time
from math import log, sqrt, pi, ceil
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
import numpy as np # noqa: E402
from scipy.special import k0, k1, k0e, exp1 # noqa: E402
from scipy.integrate import quad # noqa: E402
import principal as P # noqa: E402
import origin as O # noqa: E402
from module_law import check, section, CHECKS # noqa: E402
C1_D = 1500
CUTS = P.CUTS
J250, J10K = 0, len(CUTS) - 1
REAL_D, IMAG_D = 48000, P.IMAG_D # the real population widened (F6)
if "--d" in sys.argv: # a rehearsal: both populations cut to |D| <= this
REAL_D = IMAG_D = int(sys.argv[sys.argv.index("--d") + 1])
C1_D = min(C1_D, REAL_D)
if "--dreal" in sys.argv: # the real population alone set to |D| <= this
REAL_D = int(sys.argv[sys.argv.index("--dreal") + 1])
SEEDS = 200
NAMED = (0.0082, 0.0155) # the step's targets at |D| <= 16000 (S0), planted
DBANDS = (("|D|<=4000", 0, 4000), ("4000-8000", 4000, 8000),
("|D|>8000", 8000, 10 ** 9))
HBANDS = (("h=2", lambda h, el: h == 2), ("h=3", lambda h, el: h == 3),
("h=4", lambda h, el: h == 4), ("h>=5", lambda h, el: h >= 5),
("h>=4 el", lambda h, el: h >= 4 and el),
("h>=4 not", lambda h, el: h >= 4 and not el))
def im2(z):
"""int_z^oo K0(v) v^-2 dv: Gauss-Laguerre on k0e past 2, quad below."""
z = np.asarray(z, float)
out = np.empty_like(z)
far = z >= 2
if far.any():
u = z[far][:, None] + O.LAG_X[None, :]
out[far] = np.exp(-z[far]) * ((k0e(u) / u ** 2) @ O.LAG_W)
at2 = (np.exp(-2.0) * ((k0e(2 + O.LAG_X) / (2 + O.LAG_X) ** 2)
@ O.LAG_W))
for i in np.flatnonzero(~far):
out[i] = quad(lambda v: k0(v) / v ** 2, z[i], 2.0,
epsabs=0, epsrel=1e-13, limit=200)[0] + at2
return out
def weights2(kind, y):
z = 2 * y
if kind == "e":
return 4 * k1(z) / z + 4 * z * im2(z)
if kind == "o":
a, b = O.tails(z)
i2 = z ** 2 * k1(z) + z * k0(z) + a
return 2 * i2 / z ** 2 + 2 * z * b
return (1 + y) * np.exp(-y) / y ** 2 + np.exp(-y) - y * exp1(y)
def coeffs(G, chis, N):
"""b_n for n <= N, each character, as origin.lam1 builds them."""
b = np.ones((len(chis), N + 1), complex)
b[:, 0] = 0
for p in P.SMALL:
if p > N:
break
kr = P.kron(G.D, p)
e, pe = 1, p
while pe <= N:
idx = np.arange(pe, N + 1, pe)
idx = idx[idx % (pe * p) != 0]
b[:, idx] *= O.local(G, chis, p, kr, e)[:, None]
e, pe = e + 1, pe * p
return b
def parity(G, chis):
if G.D > 0 and G.k0 != G.e:
return ["o" if np.real(chis[i, G.k0]) < 0 else "e"
for i in range(len(chis))]
return ["e" if G.D > 0 else "c"] * len(chis)
def l2_afe(G, chis, kinds):
D = G.D
A = sqrt(D) / pi if D > 0 else sqrt(-D) / (2 * pi)
N = int(ceil((O.NA if D > 0 else 2 * O.NA) * A)) + 5
b = coeffs(G, chis, N)
n = np.arange(1, N + 1)
out = np.zeros(len(chis), complex)
for kind in set(kinds):
rows = [i for i, k in enumerate(kinds) if k == kind]
out[rows] = b[rows, 1:] @ weights2(kind, n / A)
norm = np.array([A * A * (pi / 4 if k == "o" else 1) for k in kinds])
return out / norm
def l2_euler(G, chis, broken=False):
out = np.ones(len(chis), complex)
for p in P.SMALL:
kr = P.kron(G.D, p)
q = 1.0 / (p * p)
if kr == 1:
v = chis[:, G.place(p)]
w = v if broken else np.conj(v)
out /= (1 - v * q) * (1 - w * q)
elif kr == -1:
out /= 1 - q * q
else:
out /= 1 - chis[:, G.place(p)] * q
return out
def genus_pairs(D):
"""Every unordered D = d1 d2, both fundamental, neither 1."""
out = set()
for d1 in range(-abs(D), abs(D) + 1):
if d1 in (0, 1, D) or D % d1 or not P.fundamental(d1):
continue
d2 = D // d1
if d2 != 1 and P.fundamental(d2):
out.add(tuple(sorted((d1, d2))))
return sorted(out)
def nontrivial(chis):
return [i for i in range(len(chis)) if np.max(np.abs(chis[i] - 1)) > 1e-6]
def section_controls():
section("A THE CONTROLS")
w1, n1, bad1, miss = 0.0, 0, 0, 0
w5g, w5c, n5g, n5c, brk_g, brk_c = 0.0, 0.0, 0, 0, 0.0, np.inf
P5 = P.SMALL[-1]
bar5 = 2.5 / (P5 * log(P5))
for D in list(range(5, C1_D + 1)) + list(range(-C1_D, 0)):
if not P.fundamental(D):
continue
G = P.Group(D)
if G.h < 2:
continue
chis = O.characters(G)
m, L, _ = O.lstar(G, chis)
test = [p for p in P.SMALL[1:200] if P.kron(D, p) == 1]
cols = [G.place(p) for p in test]
pairs = genus_pairs(D)
if len(pairs) != len(G.two) - 1:
miss += 1
for d1, d2 in pairs:
want = np.array([P.kron(d1, p) for p in test])
hit = [i for i in range(G.h)
if np.max(np.abs(chis[i, cols] - want)) < 1e-6]
if len(hit) != 1:
miss += 1
continue
known = ((O.l0(d1) if d1 < 0 else O.l0p(d1))
* (O.l0(d2) if d2 < 0 else O.l0p(d2)))
got = L[hit[0]]
if not (known > 0 and got.real > 0
and abs(got.imag) < 1e-10 * abs(got)):
bad1 += 1
w1 = max(w1, abs(got - known) / known)
n1 += 1
nt = nontrivial(chis)
kinds = parity(G, chis)
afe = l2_afe(G, chis[nt], [kinds[i] for i in nt])
eu = l2_euler(G, chis[nt])
br = l2_euler(G, chis[nt], broken=True)
for j, i in enumerate(nt):
rel = abs(afe[j] - eu[j]) / abs(eu[j])
relb = abs(afe[j] - br[j]) / abs(br[j])
if np.max(np.abs(chis[i].imag)) < 1e-9:
w5g, n5g, brk_g = max(w5g, rel), n5g + 1, max(brk_g, relb)
else:
w5c, n5c, brk_c = max(w5c, rel), n5c + 1, min(brk_c, relb)
check("C1w every genus character: L*(0) real, positive, the closed form",
w1 < 1e-8 and bad1 == 0 and miss == 0 and n1 > 100,
f"{n1} characters, |D| <= {C1_D}, worst relative {w1:.2e},"
f" sign or reality off at {bad1}, unmatched {miss}")
check("C5 L(2, chi) from the s = 2 kernels against the Euler product",
max(w5g, w5c) < bar5 and n5c > 100,
f"genus {n5g} worst {w5g:.2e}, complex {n5c} worst {w5c:.2e},"
f" bar {bar5:.2e}")
check("C5 break: chi(P bar) = chi(P) fails every complex character,"
" passes every genus one", brk_c > bar5 and brk_g < bar5,
f"complex least gap {brk_c:.2e}, genus worst {brk_g:.2e}")
def read_population(sign):
"""Per field (D, h, elementary, k0 present, residual T and K0 per cut)
and per character (field index, kind, genus, chi(k0), rho per cut)."""
lo, hi = (5, REAL_D + 1) if sign > 0 else (-IMAG_D, 0)
fields, chars, worst6, n6 = [], [], 0.0, 0
UNITS.clear()
UNIT_N[0] = 0
global POS_BAD, POS_N
for D in range(lo, hi):
if not P.fundamental(D):
continue
G = P.Group(D)
if G.h < 2:
continue
raw, cor = P.counts(G)
chis = O.characters(G)
m, L, _ = O.lstar(G, chis)
nt = nontrivial(chis)
Ln = L[nt]
POS_BAD += int(np.sum(~(Ln.real > 0)
| ~(np.abs(Ln.imag) <= 1e-9 * np.abs(Ln))))
POS_N += len(Ln)
kinds = parity(G, chis)
cor = np.array(cor) # (cuts, h)
tot = cor.sum(axis=1)
PIchi = cor @ chis.T # (cuts, h): PI(x, chi)
lx = np.log(np.log(np.array(CUTS, float))) + O.GAMMA
term = -np.outer(lx, m) + np.log(np.abs(L))[None, :]
rho = (PIchi - term) / tot[:, None] # (cuts, h)
hask0 = D > 0 and G.k0 != G.e
resT = np.real(rho[:, nt].sum(axis=1))
resK = (np.real((np.conj(chis[nt, G.k0])[None, :] * rho[:, nt])
.sum(axis=1)) if hask0 else None)
fi = len(fields)
fields.append((abs(D), G.h, len(G.two) == G.h, hask0, resT, resK,
tot))
local = {}
for i in nt:
chars.append((fi, kinds[i], np.max(np.abs(chis[i].imag)) < 1e-9,
np.real(chis[i, G.k0]), np.real(rho[:, i])))
gap = np.max(np.abs(chis[nt] - np.conj(chis[i])), axis=1)
key = min(i, nt[int(np.argmin(gap))])
if gap.min() > 1e-6:
key = -1 - i # no conjugate: C9u fails
if key not in local:
local[key] = UNIT_N[0]
UNIT_N[0] += 1
UNITS.append(local[key])
if abs(D) <= 3000: # C6 on the small fields
out = O.read_field(G)[0]
for j in range(len(CUTS)):
dv, E = out[j]
worst6 = max(worst6, abs(dv[G.e] - E[G.e] - resT[j]))
if hask0:
worst6 = max(worst6, abs(dv[G.k0] - E[G.k0] - resK[j]))
n6 += 1
return fields, chars, worst6, n6
POS_BAD, POS_N = 0, 0
UNITS, UNIT_N = [], [0] # per character its unit: a conjugate pair is one
def mean_se(v):
v = np.asarray(v, float)
if len(v) < 2:
return (v.mean() if len(v) else np.nan), np.nan
return v.mean(), v.std(ddof=1) / sqrt(len(v))
def fit_means(chars, j, fsel=None):
"""Kind means of rho at cut j over the characters of the selected
fields."""
acc = {}
for fi, kind, _, _, r in chars:
if fsel is None or fsel(fi):
acc.setdefault(kind, []).append(r[j])
return {k: np.mean(v) for k, v in acc.items()}
def predict(fields, chars, means_of):
"""Per field the predicted T and K0 residual from kind means;
means_of(fi) gives the means that field uses."""
pT = np.zeros(len(fields))
pK = np.zeros(len(fields))
for fi, kind, _, ck0, _ in chars:
mu = means_of(fi)[kind]
pT[fi] += mu
pK[fi] += ck0 * mu
return pT, pK
def kill(fields, chars, j, means_of, label, quiet=False):
"""The KILL FIXED test at cut j: True when it fires."""
pT, pK = predict(fields, chars, means_of)
fired, worst = False, (0.0, "")
for bn, test in HBANDS:
for cell in ("T", "K0"):
obs, pre = [], []
for fi, (a, h, el, hk, rT, rK, _) in enumerate(fields):
if not test(h, el):
continue
if cell == "T":
obs.append(rT[j]), pre.append(pT[fi])
elif hk:
obs.append(rK[j]), pre.append(pK[fi])
if len(obs) < 2:
continue
mo, se = mean_se(obs)
z = (mo - np.mean(pre)) / se
if abs(z) > abs(worst[0]):
worst = (z, f"{bn} {cell}")
fired |= abs(z) > 3
if not quiet:
print(f" {label:9s} {bn:9s} {cell:2s} observed "
f"{mo:+.4f}+-{se:.4f} predicted {np.mean(pre):+.4f}"
f" ({z:+.1f} se, {len(obs)})")
return fired, worst
def dband(a):
return next(i for i, (_, lo, hi) in enumerate(DBANDS) if lo < a <= hi)
def section_population(sign, name):
t0 = time.time()
section(f"B {name.upper()}")
fields, chars, worst6, n6 = read_population(sign)
print(f" {len(fields)} fields, {len(chars)} nontrivial characters"
f" [{time.time() - t0:.0f} s]")
check(f"C6 {name}: the characters' sum is origin.py's cell residual",
worst6 < 1e-10, f"{n6} fields with |D| <= 3000, worst {worst6:.2e}")
print(" the fields by h band x |D| band:")
for bn, test in HBANDS:
cnt = [sum(1 for a, h, el, *_ in fields
if test(h, el) and dband(a) == b) for b in range(3)]
print(f" {bn:9s} " + " ".join(f"{DBANDS[b][0]} {cnt[b]:5d}"
for b in range(3)))
print(" R1 rho per character, mean+-se (count), by kind, |D| band and"
" h band:")
for j in (J250, J10K):
print(f" cut {CUTS[j]}")
for kind in sorted({c[1] for c in chars}):
for gen in (True, False):
for b in range(3):
line = (f" {kind} {'genus ' if gen else 'complex'}"
f" {DBANDS[b][0]:10s}")
for bn, test in HBANDS[:4]:
v = [r[j] for fi, k, g, _, r in chars
if k == kind and g == gen
and dband(fields[fi][0]) == b
and test(fields[fi][1], fields[fi][2])]
if len(v) >= 2:
mo, se = mean_se(v)
line += f" {bn} {mo:+.4f}+-{se:.4f} ({len(v)})"
if "(" in line:
print(line)
print(" R2 the field residual on T against (h+ - 1)/PI(x, 1):")
for j, x in enumerate(CUTS):
y = np.array([f[4][j] for f in fields])
xs = np.array([(f[1] - 1) / f[6][j] for f in fields])
X = np.column_stack([np.ones(len(y)), xs])
beta, res, *_ = np.linalg.lstsq(X, y, rcond=None)
resid = y - X @ beta
s2 = resid @ resid / (len(y) - 2)
cov = s2 * np.linalg.inv(X.T @ X)
print(f" {x:6d} slope {beta[1]:+.4f}+-{sqrt(cov[1, 1]):.4f}"
f" intercept {beta[0]:+.5f}+-{sqrt(cov[0, 0]):.5f}")
if sign > 0:
lg = np.array([log(fields[fi][0]) for fi, k, *_ in chars if k == "o"])
ro = np.array([r[J250] for fi, k, g, ck, r in chars if k == "o"])
X = np.column_stack([np.ones(len(lg)), lg])
beta = np.linalg.lstsq(X, ro, rcond=None)[0]
resid = ro - X @ beta
cov = resid @ resid / (len(ro) - 2) * np.linalg.inv(X.T @ X)
z = beta[1] / sqrt(cov[1, 1])
print(f" R3 rho_o on log|D| at 250: slope {beta[1]:+.4f}"
f"+-{sqrt(cov[1, 1]):.4f} ({z:+.1f} se, {len(ro)} characters,"
f" taken as independent): {'HELD' if z > 3 else 'MISSED'}")
un = np.array(UNITS)
uo = np.array([k == "o" for _, k, *_ in chars])
_, first = np.unique(un[uo], return_index=True)
lgu, rou = lg[first], ro[first]
Xu = np.column_stack([np.ones(len(lgu)), lgu])
bu = np.linalg.lstsq(Xu, rou, rcond=None)[0]
ru = rou - Xu @ bu
cu = ru @ ru / (len(rou) - 2) * np.linalg.inv(Xu.T @ Xu)
print(f" R3u the same over units, a conjugate pair once (added at"
f" the audit): slope {bu[1]:+.4f}+-{sqrt(cu[1, 1]):.4f}"
f" ({bu[1] / sqrt(cu[1, 1]):+.1f} se, {len(rou)} units)")
# C7 first: the kill's control on planted rho
rng = np.random.default_rng(1660)
mu = fit_means(chars, J250)
sd = {k: np.std([r[J250] for _, kk, _, _, r in chars if kk == k],
ddof=1) for k in mu}
for plant, fn in (("fixed", lambda h, k: mu[k] + rng.normal(0, sd[k])),
("h+/4", lambda h, k: mu[k] * h / 4
+ rng.normal(0, sd[k]))):
fake = []
for fi, kind, g, ck0, r in chars:
v = r.copy()
v[J250] = fn(fields[fi][1], kind)
fake.append((fi, kind, g, ck0, v))
ff = []
for fi, f in enumerate(fields):
ff.append(list(f[:4]) + [np.zeros(len(CUTS)),
np.zeros(len(CUTS)) if f[3] else None,
f[6]])
for fi, kind, g, ck0, v in fake:
ff[fi][4][J250] += v[J250]
if ff[fi][3]:
ff[fi][5][J250] += ck0 * v[J250]
mf = fit_means(fake, J250)
fired, worst = kill(ff, fake, J250, lambda fi: mf, "", quiet=True)
want = plant == "h+/4"
check(f"C7 the kill on planted rho ({plant}) "
f"{'fires' if want else 'misses'}", fired == want,
f"worst {worst[0]:+.1f} se at {worst[1]}")
print(" KILL FIXED at 250, kind means pooled: " + ", ".join(
f"{k} {v:+.4f}" for k, v in sorted(mu.items())))
fired, worst = kill(fields, chars, J250, lambda fi: mu, "pooled")
print(f" KILL FIXED {'FIRED' if fired else 'MISSED'}: worst"
f" {worst[0]:+.1f} se at {worst[1]}")
mb = [fit_means(chars, J250, lambda fi, b=b: dband(fields[fi][0]) == b)
for b in range(3)]
print(" kind means by |D| band at 250: " + "; ".join(
f"{DBANDS[b][0]} " + ", ".join(f"{k} {v:+.4f}"
for k, v in sorted(mb[b].items()))
for b in range(3)))
fired2, worst2 = kill(fields, chars, J250,
lambda fi: mb[dband(fields[fi][0])], "by |D|")
print(f" KILL FIXED|D| {'FIRED' if fired2 else 'MISSED'}: worst"
f" {worst2[0]:+.1f} se at {worst2[1]}")
section_step(fields, chars, fired, worst)
else:
mu = fit_means(chars, J250)
print(" the imaginary arm of KILL FIXED at 250, printed, mean"
f" {mu['c']:+.4f}:")
fired, worst = kill(fields, chars, J250, lambda fi: mu, "pooled")
print(f" imaginary FIXED {'fires' if fired else 'misses'}: worst"
f" {worst[0]:+.1f} se at {worst[1]}")
def section_step(fields, chars, fired0, worst0):
"""C: the step in rho_o from h+ = 2 to h+ >= 4, priced by planted
seeds with conjugate pairs tied."""
section("C THE STEP IN RHO_O, REAL, CUT 250")
F = len(fields)
fi = np.array([c[0] for c in chars])
odd = np.array([c[1] == "o" for c in chars])
ck0 = np.array([c[3] for c in chars])
r = np.array([c[4][J250] for c in chars])
unit = np.array(UNITS)
h = np.array([f[1] for f in fields])
el = np.array([f[2] for f in fields])
hk = np.array([f[3] for f in fields])
a = np.array([f[0] for f in fields])
hc = h[fi]
sizes = np.bincount(unit)
cplx = np.array([not c[2] for c in chars])
check("C9u every complex character has its conjugate, every genus one"
" stands alone", np.all(sizes[unit] == np.where(cplx, 2, 1)),
f"{len(sizes)} units over {len(chars)} characters")
bands = [(bn, np.array([t(hh, ee) for hh, ee in zip(h, el)]))
for bn, t in HBANDS]
disj = bands[:4]
def cells(v, which):
"""z per (band, cell), pooled kind means as KILL FIXED fits."""
mo, me = v[odd].mean(), v[~odd].mean()
pre = np.where(odd, mo, me)
T = np.bincount(fi, v, minlength=F)
K = np.bincount(fi, ck0 * v, minlength=F)
pT = np.bincount(fi, pre, minlength=F)
pK = np.bincount(fi, ck0 * pre, minlength=F)
out = []
for bn, m in which:
for cell, o, p_, mm in (("T", T, pT, m), ("K0", K, pK, m & hk)):
if mm.sum() < 2:
continue
x = o[mm]
out.append((f"{bn} {cell}", (x.mean() - p_[mm].mean())
/ (x.std(ddof=1) / sqrt(len(x)))))
return out
def step(v, sel=None):
"""mean rho_o over units at h >= 4 minus at h = 2, its se."""
out = []
for m in (odd & (hc >= 4), odd & (hc == 2)):
if sel is not None:
m = m & sel
u, first = np.unique(unit[m], return_index=True)
x = v[m][first]
out.append((x.mean(), x.var(ddof=1) / len(x), len(x)))
return out[0][0] - out[1][0], sqrt(out[0][1] + out[1][1]), out
zs = cells(r, bands)
wz = max(zs, key=lambda t: abs(t[1]))
check("C9a the seeded kill reproduces KILL FIXED on the data",
abs(wz[1] - worst0[0]) < 1e-9 and (abs(wz[1]) > 3) == fired0,
f"worst {wz[1]:+.6f} at {wz[0]} against {worst0[0]:+.6f}")
t2 = r[odd & (hc == 2)].mean()
t4 = r[odd & (hc >= 4)].mean()
print(f" S0 targets, count-weighted: rho_o at h=2 {t2:+.4f}"
f" ({np.sum(odd & (hc == 2))}), at h>=4 {t4:+.4f}"
f" ({np.sum(odd & (hc >= 4))})")
mu = {True: r[odd].mean(), False: r[~odd].mean()}
sd = {True: r[odd].std(ddof=1), False: r[~odd].std(ddof=1)}
for k, nm in ((False, "even"), (True, "odd")):
line = f" S0 sd of rho {nm:4s} pooled {sd[k]:.4f}; by band"
for bn, m in disj:
sel = (odd == k) & m[fi]
if sel.sum() >= 2:
line += f" {bn} {r[sel].std(ddof=1):.4f} ({sel.sum()})"
print(line)
S, se, parts = step(r)
print(f" STEP observed {S:+.4f} +- {se:.4f} ({S / se:+.1f} se by its"
f" own se): h>=4 {parts[0][0]:+.4f} ({parts[0][2]} units),"
f" h=2 {parts[1][0]:+.4f} ({parts[1][2]})")
ab = np.array([dband(x) for x in a])[fi]
for b in range(3):
sel = ab == b
if min(np.sum(odd & sel & (hc >= 4)), np.sum(odd & sel & (hc == 2))) < 2:
continue
d, s_, pb = step(r, sel)
print(f" {DBANDS[b][0]:10s} STEP {d:+.4f} +- {s_:.4f}"
f" ({d / s_:+.1f} se; units {pb[0][2]}, {pb[1][2]});"
f" from the named step {(d - (NAMED[1] - NAMED[0])) / s_:+.1f}"
f" se (added at the audit)")
zd = cells(r, disj)
Z2 = sum(z * z for _, z in zd)
Z2T = sum(z * z for c, z in zd if c.endswith(" T"))
print(" the disjoint cells: " + ", ".join(f"{c} {z:+.2f}"
for c, z in zd))
rng = np.random.default_rng(1661)
sdk = np.where(odd, sd[True], sd[False])
base = {"FIXED": np.where(odd, mu[True], mu[False]),
"STEP": np.where(odd, np.where(hc == 2, NAMED[0], NAMED[1]),
mu[False])}
res = {}
for plant in ("FIXED", "STEP"):
fires, steps, z2s, z2ts = 0, [], [], []
for _ in range(SEEDS):
v = base[plant] + rng.normal(0, 1, len(sizes))[unit] * sdk
fires += any(abs(z) > 3 for _, z in cells(v, bands))
steps.append(step(v)[0])
zz = cells(v, disj)
z2s.append(sum(z * z for _, z in zz))
z2ts.append(sum(z * z for c, z in zz if c.endswith(" T")))
res[plant] = (fires / SEEDS, np.array(steps), np.array(z2s),
np.array(z2ts))
sdn = res["FIXED"][1].std(ddof=1)
for plant in ("FIXED", "STEP"):
fr, st, z2, z2t = res[plant]
print(f" C9 plant {plant:5s}: KILL FIXED fires {fr:.3f} of"
f" {SEEDS}; STEP z median {np.median(st) / sdn:+.2f}"
f" (sd {st.std(ddof=1):.4f}); Z2 median {np.median(z2):.1f},"
f" Z2_T {np.median(z2t):.1f}")
pw = res["STEP"][0]
print(f" prediction: KILL FIXED fires 0.10 to 0.45 under STEP:"
f" {'HELD' if 0.10 <= pw <= 0.45 else 'MISSED'} ({pw:.3f})")
check("C9 the plants: STEP z median above 1.5 under STEP, |median|"
" below 0.5 under FIXED; FIXED's STEP sd within 25% of the se;"
" KILL FIXED fires below 0.10 under FIXED",
np.median(res["STEP"][1]) / sdn > 1.5
and abs(np.median(res["FIXED"][1]) / sdn) < 0.5
and abs(sdn / se - 1) < 0.25 and res["FIXED"][0] < 0.10,
f"FIXED seeds' STEP sd {sdn:.4f} against se {se:.4f}")
for nm, obs, k in (("Z2 (T, K0; 7 cells)", Z2, 2),
("Z2_T (T; 4 cells)", Z2T, 3)):
null = res["FIXED"][k]
q95 = np.quantile(null, 0.95)
print(f" {nm}: observed {obs:.2f}, FIXED null median"
f" {np.median(null):.2f}, 95th {q95:.2f}, p"
f" {np.mean(null >= obs):.3f}; power at the 95th"
f" {np.mean(res['STEP'][k] > q95):.3f}")
z = S / sdn
verdict = "STANDS" if z > 3 else "KILLED" if z < 1 else "UNDECIDED"
print(f" KILL STEP: STEP {S:+.4f} against the FIXED seeds' sd"
f" {sdn:.4f}, {z:+.2f}: {verdict}")
d = NAMED[1] - NAMED[0]
print(f" the named step {d:+.4f}: observed {(S - d) / sdn:+.2f} sd"
f" from it")
def main():
t0 = time.time()
section_controls()
print(f" [{time.time() - t0:.0f} s]")
for sign, name in ((+1, "real, narrow"), (-1, "imaginary")):
section_population(sign, name)
check("C8 (added at the audit) L*(0, chi) real and positive at every"
" nontrivial character", POS_BAD == 0 and POS_N > 0,
f"{POS_N} characters, off at {POS_BAD}")
print(f"\n{sum(CHECKS)}/{len(CHECKS)} checks PASS, "
f"{time.time() - t0:.1f} s")
raise SystemExit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()