#!/usr/bin/env python3
"""Assemble cM_laurent.json: the data of the c_M assembly of Sec. 7 (Table 7) from the run logs of two implementations
of the eps-expansion and quadrature (retarded routing), and the leading terms with Feynman propagators on the two memory
gravitons.  Every number is copied from the named logs; nothing is recomputed except the consistency checks printed at the end.
usage: build_cM.py --logA .. --implB .. --fey1 .. --fey2 .. --feyg .. --pslq .. --moments .. --out .."""
import sys, re, json, argparse
import mpmath as mp
from common import sha256_file, dump
ap = argparse.ArgumentParser()
for k in ['logA', 'implB', 'fey1', 'fey2', 'feyg', 'pslq', 'moments', 'out']:
    ap.add_argument('--' + k, required=True)
a = ap.parse_args()
mp.mp.dps = 40
src = {}
def rd(tag, desc, path):
    src[tag] = {'description': desc, 'sha256': sha256_file(path)}
    return open(path).read()
A = rd('implementation_A_log', 'run log of implementation A (35-digit working precision, tanh-sinh quadrature of degree 5, kappa_2 by 4th-order finite differences)', a.logA)
B = json.loads(rd('implementation_B_output', 'output of implementation B (independent code, 30-digit working precision; c_M block at 25 digits)', a.implB))
F1 = rd('feynman_route_1', 'Feynman-propagator assembly from the sector integrals of implementation A', a.fey1)
F2 = rd('feynman_route_2', 'Feynman-propagator assembly from the per-domain quadratures of implementation B; closed-form candidates', a.fey2)
FG = rd('feynman_route_3', 'Feynman-propagator leading cores and eps^-4 coefficient recomputed at 33 digits', a.feyg)
PQ = rd('quadrant_moment_relations', 'integer-relation searches for the quadrant moments at two precisions', a.pslq)
MO = rd('quadrant_moment_check', 'independent double-precision quadrature of the quadrant moments', a.moments)

def grab(pat, text, g=1):
    m = re.search(pat, text, flags=re.M)
    assert m, pat
    return m.group(g).strip()
# ---- implementation A ----
secI = {}
for m in re.finditer(r'A\[(same|opp),(1|3),(1|b1|b1sq|R)\s*\] = (\S+)', A):
    secI[f'{m.group(1)},{m.group(2)},{m.group(3)}'] = m.group(4)
assert len(secI) == 16
cores = {f'j_{m.group(1)}^({m.group(2)})': m.group(3) for m in re.finditer(r'j_(1|3)\^\((0|1|2)\) = (\S+)', A)}
assert len(cores) == 6
def series(name):
    line = grab(r'^\s*' + re.escape(name) + r'[:=]?\s*=?\s*(\(.*)$', A) if False else None
def laurent(label):
    m = re.search(r'^\s*' + label + r'\s*[:=]\s*(.+)$', A, flags=re.M)
    assert m, label
    terms = re.findall(r'\(([-+0-9.e]+)\)\*eps\^(-?\d+)', m.group(1))
    return {k: v for v, k in terms}
P1, P3, c1, c3 = laurent('P1'), laurent('P3'), laurent('c1'), laurent('c3')
I1, I3, I2 = laurent('I1 ='.replace(' =', '')), laurent('I3'), laurent('I2')
def cut(d, kmin, kmax):
    return {k: v for k, v in d.items() if kmin <= int(k) <= kmax}
cM_line = grab(r'c_M\s+= (\S+)', A); cMm1 = grab(r'c_M - 1\s+= (\S+)', A); cMdig = grab(r'digits matching 1 = (\S+)', A)
I2m2 = grab(r'I2\[eps\^-2\] = (\S+)', A)
relm4 = grab(r'I2\[eps\^-4\] = \S+\s+\(rel to eps\^-2: (\S+)\)', A); relm3 = grab(r'I2\[eps\^-3\] = \S+\s+\(rel to eps\^-2: (\S+)\)', A)
pslqA = grab(r"pslq\(\[cM\]\+\[(.*?)\]\) -> (\[.*?\])", A, 2); pslqBasis = grab(r"pslq\(\[cM\]\+\[(.*?)\]\) -> (\[.*?\])", A, 1)
implA = {
 'method': 'expansion of the integrands of I_1^(M), I_3^(M) in eps through second order; two-dimensional tanh-sinh quadrature of the resulting frequency integrals (cores), split into same-sign and opposite-sign frequency sectors; exact Laurent expansion of the prefactors N_1, N_3 and of the coefficients r_1, r_2 of the integration-by-parts identity I_2 = r_1 I_3 + r_2 I_1 (Eqs. (53)-(54)); I_2^(M) from that identity',
 'working_precision_digits': int(grab(r'dps=(\d+)', A)),
 'sector_integrals': {'definition': ("A[sector,a,w]: raw two-dimensional integrals.  sector 'same': (u,v) in (0,inf)^2 with |omega_3| = u+v, integrand 2*sgn*(u v)^p K_0(u)K_0(v)K_0(u+v) * w; "
                                     "sector 'opp': (v,w) in (0,inf)^2 with u = v+w (the u > v half of one opposite-sign configuration), integrand 2*sgn*(u v)^p K_0(u)K_0(v)K_0(w) * w.  "
                                     "p = 1 for a = 1, p = 2 for a = 3; sgn = -1 except sgn(opp, a=1) = +1 (the retarded sign pattern of -omega_1^p omega_2^p summed over the two quadrants of each sector).  "
                                     "Weights w: '1'; 'b1' = -(3 ln u + 3 ln v + ln|omega_3|); 'b1sq' = b1^2; 'R' = sum_i kappa_2(|omega_i|)/K_0(|omega_i|) with kappa_2 = (1/2) d^2 K_nu/d nu^2 at nu = 0 (Eq. (56)).  "
                                     "Cores: j_a^(0) = (A[same,a,1] + 2 A[opp,a,1])/(2 pi)^2, j_a^(1) likewise with 'b1', j_a^(2) = ((A[same,a,b1sq]+2A[opp,a,b1sq])/2 + A[same,a,R] + 2A[opp,a,R] - 2 pi^2 A[same,a,1])/(2 pi)^2 (the -2 pi^2 term is the real part of the retarded phases, present in the same-sign sectors only)."),
                      'values_25_digits': secI},
 'cores': {'definition': 'I_a^(M)(eps) = N_a(eps) * sum_k j_a^(k) eps^k with N_a the prefactors of Eqs. (51), (52) (values at gamma = sqrt 2, i.e. overall powers of (gamma^2-1) stripped) and the frequency measure d omega/(2 pi) per frequency included in j', 'values_25_digits': cores,
           'exact_leading': {'j_1^(0)': '1/30', 'j_3^(0)': '-8/105', 'check_j_1^(0)_minus_1/30': grab(r'CHECK j_1\^\(0\) vs 1/30: diff = (\S+)', A)}},
 'prefactor_laurent': {'N_1': P1, 'N_3': P3, 'note': 'coefficients of eps^k, 18 significant digits printed'},
 'ibp_coefficient_laurent': {'r_2 (multiplies I_1)': c1, 'r_1 (multiplies I_3)': c3, 'note': 'Laurent expansion of the rational functions r_1, r_2 of Eq. (54), printed to 18 digits'},
 'I1_laurent': cut(I1, -1, 1), 'I3_laurent': cut(I3, -1, 1), 'I2_laurent': cut(I2, -4, -2),
 'note_orders': 'I_1, I_3 complete through eps^1 and I_2 through eps^-2 (the cores are expanded through second order); the eps^-4 and eps^-3 entries of I_2 are the numerical residuals of the pole cancellation',
 'c_M': {'value': cM_line, 'c_M_minus_1': cMm1, 'digits_matching_1': cMdig, 'I2_eps^-2': I2m2, 'definition': 'c_M = -(6 (8 pi)^4 / 5) * I_2^(M)[eps^-2]',
         'pole_cancellation_relative_to_eps^-2': {'eps^-4': relm4, 'eps^-3': relm3},
         'integer_relation_search': f'basis [c_M, {pslqBasis}] -> {pslqA} (c_M = 1, no admixture)'},
}
# ---- implementation B ----
slots = {s['name']: s for s in B['layers'][0]['slots']}
implB = {'method': 'independent implementation of the same expansion: sector decomposition of the two-frequency cores, tanh-sinh quadrature, kappa_2 by finite differences; prefactors in closed Gamma-function form',
         'working_precision_digits': 30, 'reported_digits': B['meta']['dps'],
         'I1_laurent': slots['I1M']['eps_series'], 'I3_laurent': slots['I3M']['eps_series'], 'I2_laurent': slots['I2M']['eps_series'],
         'checks': {'I1 pole vs 1/(15 (8 pi)^4), digits': round(slots['I1M']['digits']['pole_vs_closed'], 1),
                    'I1 eps^0, eps^1 vs implementation A (its 25-digit cores with exact prefactors), digits': [round(slots['I1M']['digits']['eps0_vs_oracle'], 1), round(slots['I1M']['digits']['eps1_vs_oracle'], 1)],
                    'j_3^(0) vs -8/105, digits': round(slots['I3M']['digits']['j30_vs_closed'], 1),
                    'c_M vs 1, digits': round(slots['I2M']['digits']['cM_vs_1'], 1)}}
# agreement A vs B computed here from the printed strings
agree = {}
for nm, (da, db) in {'I1': (implA['I1_laurent'], implB['I1_laurent']), 'I3': (implA['I3_laurent'], implB['I3_laurent'])}.items():
    for k in ('-1', '0', '1'):
        x, y = mp.mpf(da[k]), mp.mpf(db[k])
        agree[f'{nm}[eps^{k}]'] = float(mp.nstr(-mp.log10(abs(x/y - 1)) if x != y else mp.mpf(99), 4))
# ---- Feynman propagators on the two memory gravitons ----
def sect(text, start, stop):
    i = text.index(start); j = text.index(stop, i) if stop else len(text)
    return text[i:j]
fey = sect(F1, '--- CONSISTENT FEYNMAN', '--- raw leading moments')
rels = {m.group(1): [int(x) for x in m.group(2).split(',')] for m in re.finditer(r'^dps=24 .*basket\{x,1,pi\^2\} (\S+): \[([-\d, ]+)\]', PQ, flags=re.M)}
rels18 = {m.group(1): [int(x) for x in m.group(2).split(',')] for m in re.finditer(r'^dps=18 .*basket\{x,1,pi\^2\} (\S+): \[([-\d, ]+)\]', PQ, flags=re.M)}
assert rels == rels18 and len(rels) == 4
def closed(v):  # a x + b + c pi^2 = 0  ->  x = -(b + c pi^2)/a
    a_, b_, c_ = v
    return f"({-b_}{'+' if -c_ >= 0 else '-'}{abs(c_)} pi^2)/{a_}" if a_ > 0 else f"({b_}{'+' if c_ >= 0 else '-'}{abs(c_)} pi^2)/{-a_}"
Ssame = -mp.mpf(secI['same,1,1'])/2; Sopp = mp.mpf(secI['opp,1,1']); Tsame = -mp.mpf(secI['same,3,1'])/2; Topp = -mp.mpf(secI['opp,3,1'])
exact = {'S_same': mp.mpf(1)/3 - mp.pi**2/45, 'S_opp': mp.mpf(1)/3 + 2*mp.pi**2/45, 'T_same': 16*mp.pi**2/315 - mp.mpf(4)/9, 'T_opp': mp.mpf(4)/9 + 32*mp.pi**2/315}
num = {'S_same': Ssame, 'S_opp': Sopp, 'T_same': Tsame, 'T_opp': Topp}
qm = {}
for k in ('S_same', 'S_opp', 'T_same', 'T_opp'):
    qm[k] = {'numerical_25_digits': mp.nstr(num[k], 25), 'closed_form': {'S_same': '1/3 - pi^2/45', 'S_opp': '1/3 + 2 pi^2/45', 'T_same': '16 pi^2/315 - 4/9', 'T_opp': '4/9 + 32 pi^2/315'}[k],
             'integer_relation': f"{rels[k]} . (x, 1, pi^2) = 0 at 24 and at 18 digits", 'agreement_digits': float(mp.nstr(-mp.log10(abs(num[k]/exact[k] - 1)), 3))}
lowp = {m.group(1): m.group(2) for m in re.finditer(r'(S_same|S_opp)=([0-9.]+)', MO)}
j1F = grab(r'cand j1F = -1/90 - 1/\(3pi\^2\) = (\S+)', F2); j3F = grab(r'cand j3F =  8/315 \+ 4/\(9pi\^2\) = (\S+)', F2)
I2F4_cand = grab(r'cand I2F\[eps\^-4\] = 1/\(6144 pi\^6\) = (\S+)', F2)
feyj1 = re.search(r"^fey j1 = \['(\S+?)',", F2, flags=re.M).group(1); feyj3 = re.search(r"^fey j1 = .* j3 = \['(\S+?)',", F2, flags=re.M).group(1)
I1Fm1 = grab(r'I1\[eps\^-1\] = (\S+)\s+ratio', fey); ratio = grab(r'ratio to 1/\(15\(8pi\)\^4\) = (\S+)', fey)
I2Fm4 = grab(r'I2\[eps\^-4\] = (\S+)$', fey); I2Fm3 = grab(r'I2\[eps\^-3\] = (\S+)$', fey)
I2Fm4_G = grab(r'feybase I2\[eps\^-4\] = (\S+)', FG)
j1F_G = grab(r'j1_feybase (\S+)', FG); j3F_G = grab(r'j3_feybase (\S+)', FG)
d_routes = float(mp.nstr(-mp.log10(abs(mp.mpf(I2Fm4)/mp.mpf(I2Fm4_G) - 1)), 3))
d_exact = float(mp.nstr(-mp.log10(abs(mp.mpf(I2Fm4_G)*6144*mp.pi**6 - 1)), 3))
feyn = {
 'assignment': 'Feynman propagators on the two memory gravitons (D_31 and D_24 of arXiv:2601.16256), applied to each frequency factor as a whole: (0^+ - i s omega_j)^(a - q eps) -> (0^+ - i|omega_j|)^(a - q eps), so ln(0^+ - i|omega|) = ln|omega| - i pi/2 for either sign of omega (Eq. (57) of the paper and the text following Table 7)',
 'consequence': 'relative to the retarded assignment the contributions of the opposite-sign frequency sectors change sign at every order in eps; the leading cores and the eps^-1 pole of I_1^(M) change, and the eps^-4 pole of I_2^(M) in the identity (53) no longer cancels.  Beyond leading order the Feynman-propagator integrals are complex; only the leading (real) terms are given.',
 'quadrant_moments': {'definition': 'S_same = Int_{omega_1, omega_2 > 0} omega_1 omega_2 K_0(omega_1) K_0(omega_2) K_0(omega_1+omega_2); S_opp = the same integral of |omega_1 omega_2| over one quadrant with omega_1 omega_2 < 0; T_same, T_opp likewise with omega_1^2 omega_2^2.  Retarded: Int_{R^2}(-omega_1 omega_2) prod K_0 = -2 S_same + 2 S_opp = 2 pi^2/15 and Int(-omega_1^2 omega_2^2) prod K_0 = -2 T_same - 2 T_opp = -32 pi^2/105 (Eq. (59)); Feynman: -2 S_same - 2 S_opp and -2 T_same + 2 T_opp',
                      'values': qm, 'independent_low_precision_check': lowp},
 'leading_cores': {'j_1F^(0)': {'closed_form': '-1/90 - 1/(3 pi^2)', 'value': j1F, 'from_implementation_B_quadratures': feyj1, 'recomputed_33_digits': j1F_G},
                   'j_3F^(0)': {'closed_form': '8/315 + 4/(9 pi^2)', 'value': j3F, 'from_implementation_B_quadratures': feyj3, 'recomputed_33_digits': j3F_G},
                   'retarded_for_comparison': {'j_1^(0)': '1/30', 'j_3^(0)': '-8/105'}},
 'I1_eps^-1': {'value': I1Fm1, 'ratio_to_1/(15 (8 pi)^4)': ratio, 'closed_form_of_ratio': '30 j_1F^(0) = -1/3 - 10/pi^2'},
 'I2_eps^-4': {'value_route_1': I2Fm4, 'value_route_3': I2Fm4_G, 'routes_agree_digits': d_routes, 'closed_form': '1/(6144 pi^6)', 'closed_form_value': I2F4_cand, 'agreement_with_closed_form_digits': d_exact,
               'derivation': 'I_2[eps^-4] = r_2[eps^-3] N_1[eps^-1] j_1F^(0) + r_1[eps^-3] N_3[eps^-1] j_3F^(0) with r_2[eps^-3] = -12/5, r_1[eps^-3] = -336/5, N_1[eps^-1] = 1/(2048 pi^4), N_3[eps^-1] = 1/(131072 pi^4); the rational parts of the two terms cancel as in the retarded case and the 1/pi^2 parts of the cores leave 1/(6144 pi^6)'},
 'I2_eps^-3_real_projection': {'value_route_1': I2Fm3, 'note': 'real projection only; the full coefficient is complex'},
}
out = {'title': 'Assembly data for the constant c_M of the memory-region integrals (Sec. 7 and Table 7 of the paper): two implementations under the retarded routing, and the leading terms with Feynman propagators on the memory gravitons',
       'conventions': {'static_point': 'x = 1', 'gamma_stripping': 'all values at gamma = sqrt(2), where the overall powers of (gamma^2 - 1) equal one',
                       'measure': 'Int_omega = Int d omega/(2 pi); omega_3 = omega_1 + omega_2', 'routing': 'retarded propagators on both memory gravitons, s_1 = s_2 = +1 in Eq. (57); the all-advanced routing gives identical values, so the average of the gamma-3 prescription is automatic',
                       'identity': 'I_2^(M) = r_1(eps) I_3^(M) + r_2(eps) I_1^(M), Eqs. (53)-(54)', 'c_M': 'I_2^(M) = -5 c_M/(6 (8 pi)^4 eps^2) + O(1/eps), Eq. (16)',
                       'laurent_keys': "string k -> coefficient of eps^k (decimal strings as printed by the runs)", 'equation_numbers': 'as in the revised manuscript that this package accompanies (Sec. 7.1-7.2, Table 7)'},
       'retarded': {'implementation_A': implA, 'implementation_B': implB, 'agreement_A_vs_B_digits_from_printed_strings': agree,
                    'note_agreement': 'the comparison of the printed Laurent strings is limited by the 18 significant digits printed by implementation A; implementation B records agreement of I_1 at eps^0 and eps^1 with implementation A to 25 digits using the 25-digit cores'},
       'feynman_propagators': feyn,
       'source_sha256': src}
dump(out, a.out)
print(f"wrote {a.out.split('/')[-1]}: 16 sector integrals, 6 cores; c_M = {cM_line} (c_M-1 = {cMm1}); A-vs-B digits {agree}; "
      f"quadrant-moment closed forms agree to {[qm[k]['agreement_digits'] for k in qm]} digits; Feynman eps^-4: routes agree {d_routes} digits, = 1/(6144 pi^6) to {d_exact} digits")